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Abstract 

Most  model  reduction  techniques  employ  a  projection  framework  that  utilizes  a 
reduced-space  basis.  The  basis  is  usually  formed  as  the  span  of  a  set  of  solutions 
of  the  large-scale  system,  which  are  computed  for  selected  values  (samples)  of  in¬ 
put  parameters  and  forcing  inputs.  In  existing  model  reduction  techniques,  choosing 
where  and  how  many  samples  to  generate  has  been,  in  general,  an  ad-hoc  process.  A 
key  challenge  is  therefore  how  to  systematically  sample  the  input  space,  which  is  of 
high  dimension  for  many  applications  of  interest. 

This  thesis  proposes  and  analyzes  a  model-constrained  greedy-based  adaptive  sam¬ 
pling  approach  in  which  the  parametric  input  sampling  problem  is  formulated  as  an 
optimization  problem  that  targets  an  error  estimation  of  reduced  model  output  pre¬ 
diction.  The  method  solves  the  optimization  problem  to  find  a  locally-optimal  point 
in  parameter  space  where  the  error  estimator  is  largest,  updates  the  reduced  basis 
with  information  at  this  optimal  sample  location,  forms  a  new  reduced  model,  and  re¬ 
peats  the  process.  Therefore,  we  use  a  systematic,  adaptive  error  metric  based  on  the 
ability  of  the  reduced-order  model  to  capture  the  outputs  of  interest  in  order  to  choose 
the  snapshot  locations  that  are  locally  the  worst  case  scenarios.  The  state-of-the-art 
subspace  trust-region  interior-reflective  inexact  Newton  conjugate-gradient  optimiza¬ 
tion  solver  is  employed  to  solve  the  resulting  greedy  partial  differential  equation- 
constrained  optimization  problem,  giving  a  reduction  methodology  that  is  efficient 
for  large-scale  systems  and  scales  well  to  high- dimensional  input  spaces. 

The  model-constrained  adaptive  sampling  approach  is  applied  to  a  steady  thermal 
fin  optimal  design  problem  and  to  probabilistic  analysis  of  geometric  mistiming  in 
turbomachinery.  The  method  leads  to  reduced  models  that  accurately  represent  the 
full  large-scale  systems  over  a  wide  range  of  parameter  values  in  parametric  spaces 
up  to  dimension  21. 

Thesis  Supervisor:  Karen  E.  Willcox 

Title:  Associate  Professor  of  Aeronautics  and  Astronautics 
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Chapter  1 


Introduction 

1.1  Motivation 

Recent  years  have  seen  considerable  progress  in  solution  and  optimization  methods 
for  partial  differential  equations  (PDEs),  leading  to  advances  across  a  broad  range  of 
engineering  applications,  including  computational  fluid  dynamics  (CFD),  structural 
dynamics,  aeroelasticity,  and  large-scale  optimization,  to  name  a  few.  Improvements 
in  both  methodology  and  computing  power  have  been  substantial;  however,  a  number 
of  challenges  remain  to  be  addressed.  In  many  cases,  computational  models  for  PDEs 
lead  to  large-scale  systems  of  equations  that  are  computationally  expensive  to  solve, 
e.g.  for  applications  such  as  optimal  design  or  probabilistic  analyses. 

Model  order  reduction  is  a  powerful  tool  that  permits  the  systematic  generation 
of  cost-efficient  representations  of  large-scale  systems  that  result  from  discretization 
of  PDEs.  Several  reduction  methods  have  been  developed,  for  example,  modal  trun¬ 
cation  [2,3],  proper  orthogonal  decomposition  (POD)  [4,5],  balanced  truncation  [6], 
Krylov-subspace  methods  [7-9],  reduced  basis  methods  [10],  and  a  quasi-convex  opti¬ 
mization  approach  [11],  These  methods  have  been  applied  in  many  different  settings 
with  considerable  success,  including  controls  [12,13],  fluid  dynamics  [4,5],  structural 
dynamics  [14-19],  and  circuit  design  [20-22],  However,  a  number  of  open  issues  re¬ 
main  with  these  methods,  including  efficient  model  reduction  techniques  for  systems 
with  large  input  spaces,  and  for  nonlinear  systems. 
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Optimal  design,  optimal  control  [23,24],  probabilistic  analysis  [25]  and  inverse 
problem  applications  [26]  present  additional  challenges  for  model  reduction  methods. 
In  such  cases — where  the  physical  system  must  be  simulated  repeatedly — the  avail¬ 
ability  of  reduced  models  can  greatly  facilitate  solution  of  the  optimization  problem, 
particularly  for  large-scale  applications.  To  be  useful  for  these  applications,  the  re¬ 
duced  model  must  provide  an  accurate  representation  of  the  high-fidelity  model  over 
a  wide  range  of  parameters.  In  particular,  discretization  produces  high-dimensional 
input  spaces  when  the  input  parameters  represent  continuous  fields  (such  as  initial 
conditions,  boundary  conditions,  distributed  source  terms,  and  geometric  variability). 
Model  reduction  for  high- dimensional  input  spaces  remains  a  challenging  problem. 
Approaches  developed  for  dynamical  systems,  such  as  POD  and  Krylov-based  meth¬ 
ods,  have  been  applied  in  an  optimization  context  [21,27,28];  however,  the  number 
of  parameters  in  the  optimization  application  was  small. 


Nonlinearity  of  the  underlying  physics  of  problems  at  hand  presents  another  chal¬ 
lenge  for  projection-based  model  order  reduction.  The  difficulty  here  is  how  to  obtain 
an  efficient  reduced  model.  Even  though  the  size  of  the  reduced  model  is  much  smaller 
than  that  of  the  full  model,  in  the  nonlinear  case  it  is  not  necessarily  true  that  solving 
the  reduced  model  is  cheaper  than  solving  the  full  model.  For  example,  if  the  reduced 
matrices  depend  on  the  full  model  size  and  if  forming  (or  evaluating)  the  full  matrices 
is  the  dominant  cost,  then  solving  the  reduced  model  may  be  more  expensive  than 
solving  the  full  one  because  one  has  to  first  evaluate  the  full  matrices  before  evalu¬ 
ating  the  reduced  ones.  If,  on  the  other  hand,  the  reduced  matrices  do  not  depend 
on  the  full  model  size,  they  can  be  evaluated  once  in  the  offline  stage,  and  in  the 
online  stage,  the  cost  of  solving  the  reduced  model  is  negligible.  Therefore  a  reduced 
model  must  be  not  only  valid  for  a  range  of  parameters  but  also  independent  of  the 
full  model  size. 
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1.2  Literature  Review 


1.2.1  Projection-Based  Model  Reduction  Techniques 

Most  reduction  techniques  for  large-scale  systems  employ  a  projection  framework  that 
utilizes  reduced-space  bases.  The  key  challenge  in  projection-based  model  reduction 
is  how  to  find  a  reduced  basis  such  that  the  reduced  system  provides  an  accurate 
representation  of  the  large-scale  system  over  the  desired  range  of  inputs.  Algorithms 
such  as  optimal  Hankel  model  reduction  [29-31]  and  balanced  truncation  [6]  have 
been  used  widely  throughout  the  controls  community  to  generate  reduced  models  with 
strong  guarantees  of  quality.  These  algorithms  can  be  carried  out  in  polynomial  time; 
however,  the  computational  requirements  make  them  impractical  for  application  to 
large  systems  such  as  those  arising  from  the  discretization  of  PDEs,  for  which  system 
orders  typically  exceed  104. 

While  considerable  effort  has  been  applied  in  recent  years  towards  development  of 
algorithms  that  extend  balanced  truncation  to  large-scale  linear  time-invariant  (LTI) 
systems  [32-35],  efficient  algorithms  for  very  large  systems  remain  a  challenge.  In 
addition,  application  of  balanced  truncation  methods  to  systems  that  are  linear  time- 
varying  or  have  parametric  variation  has  been  limited  to  small  systems  [36-38].  The 
Krylov-subspace  methods  [7-9]  have  been  shown  to  be  an  alternative  efficient  model 
reduction  method  for  large-scale  LTI  systems.  Meanwhile,  the  POD  method  [4,5,39, 
40]  has  emerged  as  a  popular  alternative  for  reduction  of  very  large  dynamical  systems. 
POD  has  been  used  widely  throughout  CFD  applications  such  as  aeroelasticity  [41,42] 
and  flow  control  [27,28].  However,  both  Krylov-subspace  and  POD  methods  lack 
the  quality  guarantees  of  methods  such  as  balanced  truncation.  As  opposed  to  the 
balanced  truncation  method,  computing  a  reduced  basis  in  Krylov-subspace  and  POD 
methods  is  straightforward;  the  reduced  basis  is  formed  as  the  span  of  a  set  of  state 
solutions,  commonly  referred  to  as  snapshots.  These  snapshots  are  computed  by 
solving  the  full  system  for  selected  values  of  the  parameters  and  selected  forcing  inputs 
(possibly  selected  frequencies  if  a  Krylov-subspace  method  is  used).  The  quality  of 
the  resulting  reduced-order  model  is  very  dependent  on  the  choice  of  parameters  and 
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inputs  over  which  snapshots  are  computed.  This  is  because  the  span  of  the  snapshot 
determines  the  span  of  the  reduced  basis,  which  in  turn  determines  the  quality  of  the 
resulting  reduced-order  model.  A  key  issue  that  needs  to  be  addressed  is  therefore 
sampling;  that  is,  how  to  choose  the  parameters  and  inputs  over  which  to  compute 
the  basis. 


1.2.2  Sampling  for  Model  Reduction  in  Multi-Dimensional 
Parameter  Spaces 

In  previous  model  reduction  techniques  via  POD  or  Krylov-based  methods,  choosing 
where  and  how  many  samples  to  generate  has  been,  in  general,  an  ad-hoc  process. 
Standard  schemes  such  as  uniform  sampling  (uniform  gridding  of  the  parameter  space) 
or  random  sampling  are  not  optimal.  More  importantly,  if  the  dimension  of  the  pa¬ 
rameter  space  is  large,  uniform  sampling  will  quickly  become  too  computationally  ex¬ 
pensive  due  to  the  combinatorial  explosion  of  samples  needed  to  cover  the  parameter 
space.  Random  sampling,  on  the  other  hand,  might  fail  to  recognize  where  the  impor¬ 
tant  parameters  are  in  the  parameter  space.  One  sampling  strategy  that  compromises 
between  the  uniformity  and  the  size  of  the  sample  is  the  stratified  sampling  family 
of  which  the  popular  Latin  hypercube  sampling  (LHS)  method  is  one  example  [43]. 
The  LHS  method  is  more  efficient  than  uniform  sampling  and  often  more  accurate 
than  random  sampling.  Recently,  the  centroidal  voronoi  tessellation  (CVT)  sampling 
method  [44-46]  has  arisen  as  a  promising  method.  Initial  evaluations  [45-47]  show 
that  compared  to  LHS,  the  random  sampling  Monte  Carlo  methods,  and  Hammersley 
quasi  Monte  Carlo  sequence  methods,  on  balance  the  CVT  sampling  performs  best  at 
least  for  statistical  sampling  and  function  integration.  Nonetheless,  no  attempt  has 
compared  these  sampling  methods  in  the  context  of  efficient  and  accurate  snapshot 
generation  for  model  reduction  of  large-scale  engineering  problems. 

One  can  use  knowledge  of  the  application  at  hand  to  determine  representative 
inputs.  In  particular,  empirical  knowledge  of  the  problem  has  been  used  to  create 
the  training  parameter  set  for  the  quasi-convex  optimization  relaxation  method  [11], 
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and  to  sample  a  parameter  space  to  generate  a  POD  or  Krylov  basis  for  cases 
where  the  number  of  input  parameters  is  small,  for  example  optimal  control  ap¬ 
plications  [27,28,48,49],  aerodynamic  applications  [50,51],  physical  processes  [52], 
parametrized  design  of  interconnect  circuits  [20-22],  and  in  the  case  of  multiple  pa¬ 
rameters  describing  inhomogeneous  boundary  conditions  for  parabolic  PDEs  [53]. 
However,  this  empirical  approach  breaks  down  for  cases  with  many  parameters.  The 
recently  developed  iterative  rational  Krylov  algorithm  [54]  proposes  a  systematic 
method  for  selecting  interpolation  points  for  multipoint  rational  Krylov  approxima¬ 
tions  based  on  the  rigorous  2-norm  optimality  criterion.  This  method  has  been  ap¬ 
plied  to  reduction  of  large-scale  LTI  systems,  although  its  extension  to  parametrized 
LTI  systems  remains  an  open  question. 

Intuitively,  one  should  not  select  parameter  points  where  the  error  between  the 
full  and  the  reduced  models  is  small  since  the  reduced  model  is  already  a  good  ap¬ 
proximation  at  these  points.  As  a  result,  sampling  the  small  error  region  increases 
the  cost  but  no  new  information  is  added  to  the  snapshot  set.  Instead,  large  er¬ 
ror  parameter  points  should  be  selected  [55-57].  Recently,  a  greedy  algorithm  has 
been  proposed  to  address  the  challenge  of  sampling  a  high- dimensional  parameter 
space  to  build  a  reduced  basis  [55,56,58-62],  The  greedy  algorithm  adaptively  se¬ 
lects  snapshots  by  finding  the  location  in  a  training  parameter  set  where  an  output 
error  bound  is  maximal,  updating  the  reduced  basis  with  the  solution  at  this  sample 
location,  forming  a  new  reduced  model,  and  repeating  the  process.  The  method  has 
been  successfully  applied  to  many  applications  with  small  number  of  parameters  such 
as  fluid  dynamics  [56,63],  structure  [56,61],  heat  transfer  [55,61]  and  inverse  prob¬ 
lems  [61,62];  however,  like  the  Krylov-subspaee  and  POD  methods,  the  quality  of  the 
training  parameter  set  remains  an  open  question. 

1.2.3  Nonlinearity  Treatments  in  Model  Reduction 

Recent  efforts  towards  efficient  nonlinear  reduced  models  for  the  Euler  and  Navier- 
Stokes  equations  have  been  successful  [12,13,64],  by  exploiting  the  special  structure  of 
the  nonlinear  terms  of  the  governing  equations  so  that  reduced  models  are  completely 
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independent  of  the  full  model  size.  For  general  nonlinear  problems,  this  approach  is  no 
longer  applicable.  If  the  problems  under  consideration  are  weakly  nonlinear,  efficient 
reduced  models  could  be  obtained  by  retaining  low  order  terms  in  the  Taylor  expan¬ 
sion  of  nonlinear  terms  [65].  Another  approach  is  to  use  the  trajectory  piecewise-linear 
scheme  [66]  in  which  a  weighted  combination  of  various  linear  models  is  employed, 
hence  a  better  approximation  to  the  nonlinear  behavior  compared  with  using  a  single 
model;  however,  this  method  was  found  to  be  very  sensitive  and  not  robust  to  tuning 
parameters  for  CFD  problems  [67].  Refs.  61,62,68,69  propose  an  empirical  interpola¬ 
tion  method  in  which  the  nonlinear  terms  are  approximated  as  a  linear  combination 
of  empirical  basis  functions.  Theoretical  and  numerical  results  show  that  this  method 
is  a  promising  approach  for  a  general  nonlinear  problem. 

1.2.4  Model  Reduction  Applications 

Despite  the  fact  that  there  are  still  open  issues  remaining  to  be  addressed,  model  order 
reduction  has  been  successfully  applied  to  many  areas  of  engineering.  In  particular, 
the  reduced  basis  approach  has  been  applied  in  the  heat  transfer  context  to  design  a 
thermal  fin  to  minimize  the  material  cost  and  the  power  required  to  effect  the  desired 
cooling  [70].  In  the  context  of  active  flow  control,  the  POD  method  has  been  used  to 
control  the  wake  unsteadiness  downstream  of  flow  past  a  cylinder  [13],  and  to  control 
the  driven  velocity  in  a  driven  cavity  flow  and  the  vorticity  in  a  channel  flow  over  a 
backward- facing  step  [12].  In  order  to  address  the  fact  that  the  POD  basis  must  be 
updated  for  the  reduced  model  to  be  still  a  good  representation  of  the  full  model  when 
the  control  input  changes,  the  trust-region  POD  (TRPOD)  [71,72]  and  the  POD  for 
optimality  system  (OS-POD)  [73]  have  been  proposed  and  applied  successfully  for 
flow  control  problems.  For  applications  in  interconnect  circuit  analysis  and  MEMS, 
the  quasi-convex  optimization  approach  [11],  the  Krylov-based  methods  [20-22],  the 
trajectory  piecewise-linear  scheme  [66],  and  truncated  balanced  realization  [74]  have 
been  applied  to  reduce  full  systems  with  thousands  of  states  to  reduced  systems  with 
a  handful  number  of  states,  while  accurately  capturing  the  full  system  behavior.  In 
the  context  of  mistiming  analysis  of  bladed-disks,  both  structural  model  reduction 
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[14-19]  and  aerodynamic  model  reduction  [75]  are  able  to  accurately  represent  the 
large-scale  full  models  with  a  handful  number  of  reduced  states.  For  aeroelasticity 
applications,  both  flutter  and  limit  cycle  oscillation  predictions  have  been  successful 
with  the  POD  method  in  reducing  more  than  three  quarters  of  a  million  degrees  of 
freedom  to  a  few  dozen  degrees  of  freedom  [76,77].  Recently,  the  POD  method  has 
been  applied  to  the  CFD-based  reduced-order  aeroelastic  modeling  of  a  complete  F- 
16  fighter  configuration  [78].  The  results  show  that  the  reduced  model  aeroelastic 
predictions  are  in  very  good  agreement  with  the  full  nonlinear  model  results  and  the 
experimental  data. 


1.3  Thesis  Objectives 

The  objective  of  this  thesis  is  to  derive  a  general  framework  for  efficient  evaluation 
of  the  effects  of  parametric  inputs  (for  example,  shape  variations,  PDE  coefficients, 
initial  conditions,  etc.)  in  the  design  and  probabilistic  analysis  of  large-scale  systems, 
using  model  reduction  methods. 

In  particular,  this  thesis  aims  to: 

1.  Develop  a  systematic  technique  for  sampling  parametric  input  spaces  of  high 
dimension  in  order  to  create  a  reduced  basis. 

2.  Create  reduced  models  that  span  the  parametric  input  space  for  general  large- 
scale  systems — resulting  from  discretization  of  PDEs,  for  example — and  quan¬ 
tify  the  ability  of  the  reduced  models  to  predict  the  outputs  of  interest  over  the 
parametric  input  space. 

3.  Demonstrate  the  proposed  parametrized  model  reduction  technique  for  two 
problems:  optimal  design  of  a  steady  heat  conduction  problem,  and  probabilis¬ 
tic  analysis  of  the  effects  of  blade  shape  variations  on  unsteady  forced  response 
of  compressor  blades. 
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1.4  Thesis  Contributions 


The  main  contributions  of  the  thesis  are  to: 

1.  Develop  a  greedy-based  adaptive  model  order  reduction  approach  for  general 
parameter-dependent  problems  such  that: 

•  The  reduced  basis  takes  into  account  the  outputs  of  interest  and  the  gov¬ 
erning  equations. 

•  The  reduced  basis  spans  both  parametric  and  temporal  spaces. 

•  The  cost  of  constructing  the  reduced  basis  scales  well  to  large  dimensional 
parameter  spaces. 

2.  Propose  a  reduction  method  for  probabilistic  analysis  of  the  geometric  mistiming 
problem  in  turbomachinery 


1.5  Thesis  Outline 

In  Chapter  2,  methodologies  to  obtain  reduced  models  for  both  steady  and  unsteady 
problems  are  discussed.  In  particular,  components  of  a  general  projection-based 
model  order  reduction  method  will  be  developed.  Components  of  the  model  order 
reduction  method  that  could  make  the  approximate  solution  unbounded  are  identi¬ 
fied  and  investigated.  Then,  approaches  to  overcome  the  instability  are  discussed.  In 
Chapter  3,  a  model-constrained  greedy-based  adaptive  sampling  method  is  proposed 
for  reduction  of  large-scale  problems  that  depend  on  a  large  number  of  parameters. 
First,  the  model-constrained  adaptive  sampling  concepts,  mathematical  formulation 
and  solution  methodology  of  the  greedy  optimization  problem — which  is  one  of  the  key 
components  of  the  model-constrained  adaptive  sampling  approach — are  presented. 
An  analysis  of  the  proposed  adaptive  sampling  approach  is  then  carried  out.  Chap¬ 
ter  4  applies  the  model-constrained  adaptive  sampling  approach  and  the  projection- 
based  model  reduction  methods  on  the  steady  thermal  fin  optimal  design  problem. 
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A  number  of  numerical  simulations  are  presented  to  validate  the  theoretical  develop¬ 
ments  in  Chapter  3  and  to  compare  the  model-constrained  sampling  approach  with 
other  existing  sampling  methods.  The  application  of  the  model-constrained  adaptive 
model  reduction  method — which  is  the  combination  of  the  mo  del- constrained  adaptive 
sampling  approach  proposed  in  Chapter  3  and  the  projection-based  model  reduction 
techniques  discussed  in  Chapter  2 — in  the  context  of  optimal  design  is  then  discussed. 
To  prepare  for  the  numerical  results  in  Chapter  6,  Chapter  5  presents  an  unsteady 
linearized  CFD  model  based  on  the  Discontinuous  Galerkin  finite  element  method.  A 
linearized  CFD  model  for  incorporating  geometric  variation  effects  into  the  unsteady 
simulation  is  then  developed.  Numerical  results  to  validate  the  linearized  CFD  model 
are  also  discussed.  Chapter  6  begins  with  a  numerical  result  to  demonstrate  the  sta¬ 
bility  and  the  convergence  of  the  Petrov-Galerkin-projection-based  model  reduction 
discussed  in  Chapter  2.  The  application  of  the  model-constrained  adaptive  model 
reduction  approach  to  find  reduced  models  for  an  unsteady  problem  in  turboma¬ 
chinery  application  is  then  presented.  In  particular,  the  application  of  the  resulting 
reduced  models  in  the  context  of  probabilistic  analysis  with  a  large  number  of  geo¬ 
metric  variability  parameters  is  investigated.  Finally,  Chapter  7  concludes  the  thesis 
with  recommendations  for  extensions  and  future  work. 
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Chapter  2 


Projection-Based  Model  Order 
Reduction 


In  this  chapter  we  discuss  methodologies  to  obtain  reduced  models  for  both  steady 
and  unsteady  problems.  In  particular,  the  chapter  begins  by  developing  components 
of  a  general  projection-based  model  order  reduction  method.  We  first  investigate  why 
a  projection-based  reduced  model  could  be  unstable.  That  is,  we  identify  components 
of  the  model  order  reduction  method  that  could  make  the  approximate  solution  un¬ 
bounded.  Then  we  discuss  approaches  to  overcome  the  instability.  One  of  the  main 
ingredients  is  the  minimum- residual  (least-squares)  projection  that  is  widely  used 
in  linear  algebra  [79],  in  the  finite  element  community  [80-82],  and  recently  in  the 
reduced  basis  context  [83-85]. 


2.1  General  Projection-Based  Model  Order  Reduc¬ 
tion 

Most  large-scale  model  reduction  frameworks  are  based  on  a  projection  approach, 
which  can  be  described  in  general  terms  as  follows.  Given  a  large-scale  dynamical 
problem 
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Find  x  G  Mn  such  that 


R(x,  x,  z,  u,  t)  =  0,  x(0)  =  x°,  y  =  P(x,z,u,i),  (2.1) 

where  x  =  x(z,  u,  t)  G  Rn  is  the  full  solution  (state)  vector,  z  G  Rd  is  the  vector 
containing  the  parameters  of  interest,  t  denotes  time,  x°  is  the  specified  initial  state, 
u  =  u (t)  G  Rp  is  some  time-dependent  input  vector,  R  is  some  discrete  operator 
(residual  operator  resulting  from  a  numerical  discretization  of  a  set  of  PDEs,  for 
example),  and  y  G  R9  is  a  vector  containing  q  outputs  of  interest  computed  by 
some  output  operator  V.  The  dynamical  system  (2.1)  could  be,  for  example,  the 
finite  element  discretization  of  a  set  of  PDEs,  or  governing  equations  in  molecular 
dynamics  simulation,  or  a  dynamical  system  from  circuit  simulation,  etc. 

We  will  develop  projection-based  model  order  reduction  techniques  using  linear 
algebra  tools  [79,86,87].  The  functional  analysis  point  of  view,  i.e.  using  variational 
tools,  can  be  found  in  Refs.  56,61,62,84.  Since  the  size  of  the  system,  n,  is  typically 
very  large,  e.g.  n  >  105,  in  the  context  of  design  and  optimization — in  which  the 
large-scale  problem  (2.1)  needs  to  be  solved  repeatedly  for  many  different  design 
parameters  z — it  is  too  computationally  expensive  to  use  the  original  full  problem 
(2.1).  Instead,  we  seek  an  approximate  solution  x  within  a  reduced  space  Vm.  That 
is,  x  G  Vm  C  Rn  where  Vm  is  defined  to  be  the  span  of  some  m  independent  vectors 
of  Rn,  and  is  called  the  trial  space.  We  also  introduce  the  test  space  Wm  C  Rn.  Let 
<E»,  G  Rnxm  be  a  basis  of  Vm  and  Wm,  respectively,  and  assume  VI/  /  $  =  I,  where 
I  is  the  identity  matrix.  For  example,  $  contains  as  columns  the  basis  vectors  0j, 
i.e.,  $  =  [0i  02  •  •  •  0m];  and  is  assumed  to  span  the  full  solution  space  over  a  range 
of  interest  of  parameters  z  and  unsteady  inputs  u,  i.e.  x(z,  u,  t)  G  span{<L},  Vz  G 
5z,Vu  G  5U  where  Sz  and  Su  are  some  subspaces  of  interest  of  the  parametric  and 
forcing  spaces,  respectively.  The  approximate  solution  x  can  be  then  expressed  as 

x  =  <3>xr,  (2.2) 

and  the  reduced  model  can  be  written  in  matrix  form 
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Find  xr  e  Mm  such  that 


^TR(#xr,  $xr,  z,  u,  t)  —  0,  xr(0)  =  \E,i  x°,  yr  =  V(&xr,  z,  u,  t),  (2.3) 

where  xr  =  xr(z,u,  t)  E  Hm  is  the  vector  of  the  coordinates  of  the  approximate 
solution  x  in  the  reduced  basis  <h,  and  yr  E  1R9  is  the  approximate  output.  We  have 
used  the  reduced  transform  (2.2)  and  the  assumption  th1  $  =  I  to  obtain  the  initial 
condition  for  the  reduced  state  xr(0).  If  the  test  space  is  the  same  as  the  trial  space, 
i.e.  th  =  3>,  the  reduced  system  (2.3)  is  obtained  via  a  Galerkin  projection.  If  the 
test  space  is  different  from  the  trial  space,  the  reduced  system  (2.3)  is  obtained  via  a 
Petrov- Galerkin  projection. 

Figure  2-1  summarizes  the  above  projection-based  reduction  method.  One  starts 
with  a  dynamical  system  in  the  full  space.  The  first  step  is  to  identify  the  reduced 
basis  pair  $,  \F,  and  hence  the  reduced  spaces  Vrn  and  Wm.  The  second  step  is 
to  employ  the  reduced  transform  (2.2)  and  to  perform  the  projection  to  obtain  the 
reduced  dynamical  system.  Once  the  reduced  state  xr  is  found,  the  approximate 
solution  is  then  reconstructed  using  the  reduced  transform  (2.2). 

Reduced 


Figure  2-1:  A  general  projection-based  model  order  reduction. 

One  of  the  important  tasks  of  a  projection-based  model  reduction  technique  is 
therefore  to  find  a  reduced  basis  pair  U/  and  $  so  that  the  reduced  system  (2.3) 
provides  an  accurate  representation  of  the  large-scale  system  (2.1)  over  the  desired 
range  of  inputs  and  (possibly)  parameters. 


39 


Next  let  us  apply  the  general  projection  framework  to  find  a  reduced  model  of  the 
general  parametrized  LTI  dynamical  system 

E(z)x  =  A(z)x  +  B(z)u,  y  =  C(z)x,  (2.4) 

with  initial  condition 

x(0)  =  x°,  (2.5) 

where  the  matrices  E(z),A(z)  G  Rnxn,  B(z)  G  Rnxp,  and  C(z)  G  Hqxn  in  (2.4) 
may  depend  (possibly  nonlinearly)  on  a  set  of  parameters  z;  the  parameters  z  could 
be  coefficients  of  the  PDEs,  for  example  heat  conductivities,  or  shape  parameters. 
The  input  vector  u(f)  G  Rp  could  describe  prescribed  unsteady  boundary  motion. 
Systems  of  the  form  (2.4)  that  result  from  spatial  discretization  of  a  set  of  PDEs  can 
be  found  in  Chapter  5.  In  this  case,  the  dimension  of  the  system,  n,  is  very  large  and 
the  matrices  E(z),  A(z),  B(z)  and  C(z)  result  from  the  chosen  spatial  discretization 
method. 

Now  denoting  the  residual  as 

R(<E*x,r,  3?xr,  z,  u,  t )  =  E(z)<hx,r  —  A(z)3?xr  —  B(z)u,  (2-6) 

we  apply  the  above  general  projection-based  model  order  reduction  technique.  This 
yields  the  parametrized  LTI  reduced-order  model  with  state  xr(t)  and  output  yr(t) 


E,.(z)xr 

=  A,r(z)x,r  +  Br(z)u, 

(2.7) 

y-r 

C,r(z)xr, 

(2.8) 

xr(0) 

=  *Tx°, 

(2.9) 

where  Er(z)  =  $JE(z)$,  Ar(z) 

=  TfTA(z)$,  Br(z)  =  ^TB(z), 

C,(z)  =  C(z)#. 

Note  that  we  have  performed  the  projection  without  discretizing  the  time  derivative 
terms,  but  one  could  also  discretize  the  time  derivative  terms  first,  then  perform  the 
projection  on  the  fully  discrete  system. 
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If,  on  the  other  hand,  one  starts  with  a  full  steady  system 


A(z)x  —  B(z),  y  —  C(z)x, 


(2.10) 


and  defines  the  residual  as 


R(<hx.r,z)  =  B(z)  —  A(z)<3>x.r,  (2-11) 

the  above  projection-based  model  order  reduction  technique  will  yield  the  reduced 
system  of  the  form 

Ar(z)xr  =  Br(z),  yr  =  Cr(z)xr,  (2.12) 

where  again  Ar(z)  =  ^,IA(z)$.  Br(z)  =  \l>TB(z),  Cr(z)  =  C(z)<3>.  Systems  of  the 
form  (2.10)  that  result  from  spatial  discretization  of  a  set  of  PDEs  can  be  found  in 
Chapter  4. 

In  the  following  we  discuss  methodologies  to  construct  the  reduced  basis  th  to 
obtain  a  guaranteed  stable  reduced  model  for  both  steady  and  unsteady  problems 
that  are  linear  in  state  vector  as  in  (2.4)-(2.5)  and  (2.10).  Extensions  to  problems 
that  are  nonlinear  in  state  vector  are  not  addressed  in  this  thesis.  Construction  of 
the  basis  $  will  be  discussed  in  Chapter  3. 


2.2  Construction  of  the  Reduced  Test  Basis  \I/  for 
Steady  Problems 

For  steady  problems,  we  first  denote  e  =  x— x,  the  state  error  between  the  full  solution 
and  the  approximation.  The  A1  A— norm  is  defined  as  ||  v||  ata  =  v1  A1  Av,  Vv  e  IRC 
From  the  residual-error  relation  R(<hxr,z)  =  Ae,  where  R  is  defined  in  (2.11),  it  is 
easy  to  see  that  the  following  identity  is  true 

||R||2  =  ||e||ATA.  (2.13) 
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Now  if  the  Galerkin  projection,  i.e.  U/  =  T>,  is  employed  to  find  the  reduced 
state  xr,  it  provides  no  guarantee  on  the  stability  of  the  reduced  model  if  A  is  not 
symmetric  positive  definite.  (If  A  is,  however,  symmetric  positive  definite  one  could 
prove  that  the  Galerkin  projection  yields  an  optimal  reduced  model  in  the  A— norm, 
that  is,  the  state  error  is  minimized  in  the  A— norm.)  This  is  because  the  reduced 
equation  (2.12)  is  obtained  by  enforcing  the  residual  to  be  orthogonal  to  the  reduced 
space,  and  hence  the  residual  could  be  mathematically  arbitrarily  large  while  being 
orthogonal  to  the  reduced  space. 

The  above  discussion  suggests  that  the  approximate  solution  should  minimize  the 
residual  in  (2.11).  In  other  words,  we  find  the  reduced  state  xr  to  minimize  the 
residual.  This  is  a  form  of  the  minimum-residual  statement  in  linear  algebra  and 
the  finite  element  contexts  [79-81].  Recently,  the  minimum- residual  approach  has 
been  successfully  used  in  the  reduced-basis  context  [84,85,88].  Mathematically,  the 
minimum-residual  statement  can  be  expressed  in  terms  of  the  following  least  square 
minimization 

xr  =  arg  min  ||R(«I?xr,  z) || |  =  ||B  —  A<3>xr||2  (2-14) 

xrelRm 

whose  optimality  condition,  which  is  the  reduced  model,  is  given  by 

(A$)T(A$)xr  =  (A$)tB.  (2.15) 

Next,  if  we  choose  the  reduced  test  basis  to  be 

T'  =  A3>,  (2.16) 

then  the  minimum-residual  approach  is  equivalent  to  a  Petrov-Galerkin  projection 
with  the  test  space  given  in  (2.16).  This  particular  Petrov-Galerkin  projection  is 
equivalent  to  Galerkin  projection  on  the  normal  equation.  Moreover,  the  fact  in 
(2.13)  that  the  residual  2-norm  is  exactly  the  state  error  in  the  A7  A— norm  makes 
this  particular  Petrov-Galerkin  projection  the  best  reduction  approach  in  the  sense 
that  the  resulting  reduced  model  minimizes  the  state  error  in  the  A1  A— norm.  It 
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should  be  pointed  out  that  the  minimum-residual  approach  yields  the  test  reduced 
basis  and  the  reduced  model  at  the  same  time.  We  are  now  in  position  to  prove  the 
reduced  model  (2.15)  obtained  from  the  minimum-residual  statement,  or  equivalently 
from  the  Petrov-Galerkin  projection  (2.16),  is  guaranteed  to  be  stable. 

Theorem  2.1  Assume  A  has  a  bounded  condition  number,  then  the  reduced  model 
(2.15)  is  stable  in  the  sense  that  the  state  error  is  bounded.  In  particular,  the  following 
bound  holds 

II e|| 2  <  — x —  l|R(^*xr,  z) || 2  (2.17) 

(J 

mm 

where  cWn  is  the  smallest  singular  value  of  A.  the  reduced  state  xr  is  computed  from 
equation  (2.15),  and  the  residual  R  is  defined  in  (2.11). 

Proof:  Making  use  of  the  inequality  for  compatible  matrix  and  vector  norms  on 
the  error-residual  relation  e  =  A  *R  and  using  the  definition  of  the  singular  val¬ 
ues  of  a  matrix  yield  the  bound  (2.17).  Since  the  reduced  state  found  from  equa¬ 
tion  (2.15)  minimizes  the  residual,  the  residual  is  finite  (because  ||R(xr,z)||2  = 
min_  g]p>™  ||R(<f>x,r,  z)||2  <  ||R(<E»0,  z)||2  =  ||B(z) ||2).  In  addition,  l/cr^in  is  finite 
due  to  the  bounded  conditioned  number  assumption  on  A.  We  therefore  conclude 
that  the  state  error  in  the  2-norm  is  bounded.  □ 

A  priori  convergence  results  are  standard  in  the  context  of  the  finite  element 
method  [2,89],  and  recently  in  the  reduced  basis  approach  context  [56,58-62],  In  this 
light,  we  present  an  a  priori  convergence  result  for  the  above  steady  reduced  model. 

Theorem  2.2  As  the  reduced  basis  $  is  enriched,  i.e.  more  basis  vectors  are  added, 
the  approximate  solution,  hence  the  reduced  model,  is  improved  in  the  sense  that 
the  state  error  in  the  A 1  A— norm  is  a  non-increasing  function  of  the  number  of 
reduced  basis  vectors.  In  particular,  there  exists  m  <  n  at  which  the  state  error  in  the 
ATA  —  norm  is  strictly  monotone  decreasing  (a  linear  algebra  version  of  the  a  priori 
convergence  result  in  [56,  58-62]). 

Proof:  Assume  $  G  Rnxm  is  the  current  reduced  basis.  Now  adding  a  new  basis 
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vector  <pm+ 1,  the  new  reduced  basis  <f>  is  given  as 

$  =  [$,  (j)m+ 1]  G  Rnx(m+1).  (2.18) 

The  minimum-residual  statement  of  the  new  reduced  model  is  given  as 

x*  =  arg  min  ||B-A[#,0m+i]xr||2.  (2.19) 

x,..-  Rm+1 

Comparing  the  minimization  statements  (2.14)  and  (2.19),  and  using  the  fact  that 
span{$>}  C  span{<&}  and  Rm  C  Rm+1,  it  can  be  seen  that  xr  is  a  special  case  of 
xr  whose  last  element  is  zero.  That  is,  in  (2.19),  the  residual  norm  is  minimized 
in  a  larger  space,  and  hence  the  residual  should  be  no  larger  than  that  of  (2.14). 
Equivalently,  the  state  error  in  the  A7  A— norm  is  no  larger  when  the  reduced  space 
is  richer.  Since  the  approximate  solution  is  exact  if  m  =  n,  there  exists  m  <  n  such 
that  when  more  basis  vectors  are  added  the  state  error  is  smaller.  □ 

Note  that  from  the  residual-error  relation  R  =  Ae  and  the  output  error  relation 
y  —  yr  =  Ce,  we  conclude  that  as  the  reduced  basis  is  enriched,  the  residual  and  the 
output  error  are  also  improved.  This  a  priori  convergence  result  is  important  for  the 
theoretical  development  of  our  adaptive  sampling  method  discussed  later  in  Chapter 
3. 

2.3  Construction  of  the  Reduced  Test  Basis  \I/  for 
Unsteady  Problems 

For  the  unsteady  problem,  we  first  review  mathematically  why  a  reduced  model  could 
be  unstable  even  if  the  full  model  is  stable.  The  details  of  the  origin  of  the  instability 
can  be  found  in  [86,90]  and  references  therein.  To  begin,  assume  E  =  I  and  denote 
the  numerical  range  of  the  matrix  A  as 

T( A)  =  {xrAx  :  || x|| 2  =  1}.  (2.20) 
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If  A  is  normal,  i.e.  AA1  =  A1  A,  the  numerical  range  is  equal  to  the  convex  hull 
of  the  spectrum  of  A.  As  a  consequence,  the  reduced  matrix  $7A$  can  never  be 
unstable,  e.g.  the  reduced  solution  xr  does  not  increase  exponentially  in  time,  if  A 
is  stable.  However,  if  A  is  not  normal,  the  numerical  range  can  extend  into  the  right 
half  plane  (e.g.  the  reduced  matrix  has  some  eigenvalues  with  positive  real  parts) 
depending  on  the  reduced  basis  <E>,  and  hence  the  reduced  matrix  could  be  unstable. 
Unfortunately,  in  most  CFD  applications,  the  matrix  A  is  not  normal  and  that  is  the 
reason  why  the  Galerkin  projection  method  can  yield  an  unstable  reduced  model. 

It  is  well  known  that  there  two  approaches  in  optimization.  That  is,  one  can 
use  either  the  differentiate-then- discretize  or  discretize-then-differentiate  approaches. 
Each  of  these  methods  has  its  own  advantages  and  disadvantages;  a  detailed  discussion 
of  these  two  approaches  in  the  optimal  control  context  can  be  found  in  [91].  In 
the  following,  we  discuss  the  application  of  the  discretize-then-differentiate  approach 
together  with  the  minimum-residual  statement  in  Section  2.2  to  find  a  test  reduced 
basis  \I/  and  at  the  same  time  a  reduced  model  for  time-dependent  problems. 

In  the  discretize-then-differentiate  approach,  one  first  discretizes  the  time  deriva¬ 
tive  of  the  full  model.  A  suitable  optimization  problem  is  then  constructed  and 
differentiated  to  find  the  optimality  condition.  The  optimality  condition  is  solved  to 
find  an  approximate  optimizer. 

Now,  without  loss  of  generality,  assume  that  the  Backward-Euler  method  is  used 
to  discretize  the  time-dependent  terms  of  equation  (2.4),  the  full  model  is  now  given 
as 

(E  -  AtA)x.k  =  Ex1'-1  +  AfBufc,  k  =  1, . . . ,  nt,  (2.21) 

where  k  denotes  the  time  level,  At  is  the  time  step,  and  nt  is  the  number  of  time 
steps.  Next,  denote  the  residual  at  time  step  k  as 

Rfc  =  (E  -  AtA)±k  -  Exfc_1  -  AtBuk,  (2.22) 

where  xfc  =  t&x]?  is  the  approximation  at  the  kth  time  step,  and  x.k  is  the  corresponding 
reduced  state  vector. 
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Similar  to  the  minimum-residual  statement  in  Section  2.2,  all  the  residuals  are 


minimized  simultaneously  as 


X 


1 

r ; 


nt  nt 

min  E  HR‘H2  =  E  II  <E  -  -  Efep1 

•x"  k= 1  k= 1 


AtHuk\\l. 


(2.23) 


This  is  equivalent  to  stacking  up  the  residuals  at  all  time  steps  into  a  single  long 
residual  vector.  Symbolically,  the  residual  equations  can  be  written  as 


R1 

-AfBu1  -  E$x° 

Q 

$ 

Xr 

R2 

= 

-AfBu2 

— 

E  Q 

$ 

x2 

R' 

— AtBunt 

E  Q 

$ 

R 


(2.24) 

where  Q  =  —  (E  —  At  A),  and  blank  entries  in  matrices  A  and  $  mean  zero  blocks. 
The  minimum-residual  problem  (2.23)  is  equivalent  to  the  following  minimum-residual 
problem 


min  || R|| |  =  ||B  —  A3>xr||2,  (2.25) 


where  R,  B,  A  and  $  are  defined  in  (2.24).  Next,  setting  the  first  derivative  of  the 
residual  norm  || R|||  with  respect  to  x(r  to  zero,  we  obtain  the  reduced  equations 


[(E-  AtA)$]rRfc-  (E<S>)TRk+1  =  0,  k  —  1, ...,  nt  —  1, 
[(E-  AtA)$]TRnt  =  0, 


(2.26) 

(2.27) 


which  can  be  viewed  as  a  weighted  Petrov- Galerkin  projection  in  which  two  different 
test  reduced  spaces  (E  —  AfA)«h  and  E<f>  are  used  for  the  residual  at  two  successive 
time  steps.  In  matrix  form,  the  reduced  equations  have  the  following  symmetric  block 
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tri-  diagonal  form 


H  Dt  0  . 

Xr 

F1 

D  H  Dt  0  . 

p2 

o  . 

...  0  D  H  Dt  0  ... 

= 

jpfc 

.  0  D  H  Dr 

-v-nt— 1 

.  0  D  H 

Ynt 

(2.28) 


where 


H  =  [(E  -  AfA)$]T  [(E  -  AfA)<E>]  +  $TETE$,  (2.29) 

H  =  [(E-  AtA)$]T[(E-  AfA)$],  (2.30) 

D  =  -[(E-  AtA)$]T(E$),  (2.31) 

F1  =  [(E-  AfA)$]T(E$x°  +  AtBu1)  -  (E#)T(AfBu2),  (2.32) 

Fk  =  [(E  -  Af  A)$]T  (AfBufc)  -  (E<f>)T(AfBufc+1),  k  =  2, ...,  nt  -  1,  (2.33) 

Fnt  =  [(E-  (AtBunt).  (2.34) 


It  can  be  seen  that  the  residual  equation  (2.24)  and  the  minimum-residual  state¬ 
ment  (2.25)  have  exactly  the  same  form  as  the  steady  residual  equation  (2.11)  and 
the  minimum- residual  statement  (2.14),  respectively.  As  a  result,  all  the  results  in 
Section  2.2  hold  for  reduced  models  of  unsteady  problems  using  the  method  in  this 
section  as  well.  In  particular,  the  stability  of  the  form  (2.17)  and  a  priori  convergence 
result  as  in  Theorem  2.2  hold. 

However,  even  though  the  system  (2.28)  is  sparse,  its  size,  m  x  nt,  can  be  large, 
because  even  if  the  number  of  reduced  states  m  is  moderate,  the  number  of  time 
steps  nt  could  be  very  large.  In  that  case,  an  efficient  algorithm  must  be  derived 
to  solve  this  symmetric  block  tri-diagonal  system  to  ford  the  reduced  state  vectors. 
Therefore,  while  this  method  is  of  theoretical  interest  (because  it  provides  an  a  priori 
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convergence  result  which  is  useful  for  the  theoretical  developments  of  our  adaptive 
sampling  approach  in  Chapter  3),  in  practice  we  minimize  the  residuals  sequentially. 
That  is,  the  residual  at  the  first  time  step  is  minimized  to  determine  the  reduced  state 
at  the  first  time  step.  Once  the  reduced  state  at  the  first  time  step  is  computed,  the 
residual  at  the  second  time  step  is  minimized  to  determine  the  reduced  state  at  the 
second  time  step.  This  process  is  repeated  until  the  reduced  state  at  the  final  time 
step  is  computed.  In  particular,  minimizing  the  residual  at  the  kth  time  step  yields 
the  corresponding  reduced  equation  at  that  time  step  as 

[(E  -  At A)<E>]r  Rfc  =  0  (2.35) 

which  can  be  seen  to  be  obtained  via  a  Petrov- Galerkin  projection  with  the  test 
reduced  basis  as  \1/  =  (E  —  AfA)<E>.  The  main  advantage  of  this  approach  is  that 
the  reduced  equations  at  each  time  step  are  decoupled,  and  hence  are  easy  to  solve. 
We  do  not  yet  have  a  priori  convergence  result  as  in  Theorem  2.2,  though  numerical 
results,  as  we  will  present  in  Section  6.1,  justify  this  convergence. 
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Chapter  3 


Model-Constrained  Greedy-Based 
Adaptive  Sampling  Approaches  for 
Parameter-Dependent  Problems 


We  have  discussed  methodologies  to  compute  a  good  test  reduced  basis  in  Chapter 
2.  In  this  chapter  we  propose  approaches  to  find  a  good  trial  reduced  basis.  In  par¬ 
ticular,  we  develop  a  model-constrained  greedy-based  adaptive  sampling  method  for 
model  reduction  of  large-scale  problems  that  depend  on  a  large  number  of  param¬ 
eters.  First,  the  model-constrained  adaptive  sampling  concepts,  and  corresponding 
mathematical  formulations  are  presented.  A  solution  methodology  for  the  greedy 
optimization  problem,  which  is  one  of  the  key  components  of  the  model-constrained 
adaptive  sampling  approach,  is  also  discussed.  An  analysis  of  the  proposed  adaptive 
sampling  approach  is  then  carried  out.  Finally,  the  greedy  optimization  problem  and 
its  first  order  optimality  conditions  are  derived  for  steady  and  unsteady  problems 
that  are  linear  in  the  state  vector. 
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3.1  Model-Constrained  Greedy-Based  Sampling  Ap¬ 


proach 

Recall  that  our  model  reduction  task  is  one  of  determining  an  appropriate  reduced 
basis  that  spans  both  the  parametric  input  space  z  and  the  space  of  unsteady  inputs 
u(f).  In  the  case  of  the  dynamical  system  (2.1)  with  no  dependence  on  parameters  z, 
a  number  of  model  reduction  techniques  can  be  used,  such  as  Krylov-based  methods 
and  POD.  In  these  methods,  the  reduced  basis  is  formed  as  the  span  of  a  set  of 
state  solutions,  commonly  referred  to  as  snapshots.  These  snapshots  are  computed 
by  solving  the  full  system  for  selected  values  of  the  parameters  and  selected  forcing 
inputs  (possibly  selected  frequencies  if  a  Krylov-subspace  method  is  used).  In  order 
to  extend  these  techniques  to  the  general  case  where  the  system  matrices  depend  on 
the  parameters  z,  we  require  a  systematic  method  of  sampling  the  parametric  input 
space,  and  the  forcing  input  space  as  well. 

3.1.1  Greedy  Adaptive  Sampling  Concept 

A  recently  proposed  approach  to  address  the  challenge  of  sampling  a  high-dimensional 
parameter  space  to  build  a  reduced  basis  is  the  greedy  algorithm  [56,58,62,63]. 
The  greedy  algorithm  adaptively  selects  snapshots  by  finding  the  location  in  a  pre¬ 
determined  discrete  parameter  set  (training  parameter  set)  where  an  output  error 
bound,  i.e.  an  upper  bound  of  the  error  between  the  full  and  the  reduced  outputs, 
is  maximal,  updating  the  basis  with  information  gathered  from  this  sample  location, 
forming  a  new  reduced  model,  and  repeating  the  process. 

Here,  we  formulate  the  greedy  approach  as  an  optimization  problem  that  targets 
an  error  estimation  (which  could  be  an  output  error  indicator,  an  output  error  bound, 
or  the  true  output  error)  of  reduced  model  output  prediction.  The  optimization  prob¬ 
lem  is  defined  by  introducing  as  constraints  the  systems  of  equations  representing  the 
reduced  model  (and  possibly  the  full  model  if  the  true  output  error  is  targeted).  The 
optimization  formulation  treats  the  parameter  space  as  continuous;  that  is,  we  do  not 
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require  a  priori  selection  of  a  discrete  training  parameter  set.  As  a  result,  instead 
of  finding  globally-optimal  parameter  points  in  the  training  set,  our  method  seeks  to 
find  locally-optimal  parameter  points  in  the  continuous  parameter  space.  Further, 
since  any  error  estimation  can  be  used  as  our  selection  criteria,  our  approach  is  appli¬ 
cable  in  cases  for  which  output  error  bounds  are  unavailable.  We  use  state-of-the-art 
optimization  techniques  to  solve  the  resulting  greedy  PDE-constrained  optimization 
problem. 


3.1.2  Greedy  Optimization  Problem 


In  each  cycle  of  the  greedy  algorithm,  the  key  step  is  to  determine  the  location  in 
parameter  space  where  the  error  in  the  reduced  model  is  maximal.  For  the  sake 
of  clarity,  we  first  discuss  sampling  methods  in  the  parameter  space,  and  then  ad¬ 
dress  sampling  approaches  for  the  unsteady  forcing  input  space.  We  define  the  cost 
functional  as  a  function  of  the  output  error  norm 

Q  (x,  xr ,  z)  =  ~  ||  y  yv  ||  q,  (3.1) 


where  the  appropriate  definition  of  the  output  error  norm  ||  ■  ||o  for  steady  and  un¬ 
steady  problems  will  be  discussed  later.  Given  a  current  basis  <F,  we  find  the  location 
in  parameter  space  of  maximum  output  error  by  solving  the  optimization  problem 


max  Q 

X,Xr,Z 


2  y 


||  2 
llo 


(3.2) 
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subject  to 


R(x,  x,  z,  u(f),  t) 

=  0, 

(3.3) 

x(0) 

(3.4) 

y 

=  P(x,Z,U  (t),t), 

(3.5) 

tl/7  R(<E»xr,  3>xr,  z,  u(f),  t ) 

=  0, 

(3.6) 

xr(0) 

=  *Tx°, 

(3.7) 

Yr 

=  V(&Xr,Z,u(t),t), 

(3.8) 

Z min  — 

z  A  zmaX) 

(3.9) 

where  z min  and  zmax  are  respectively  lower  and  upper  bounds  on  the  parameter  vector 
z.  Since  the  true  output  error — the  error  between  the  full  and  the  reduced  outputs — 
has  been  used  as  the  cost  functional,  both  the  full  model  (3.3)-(3.5)  and  the  reduced 
model  (3.6)-(3.8)  are  constraints  of  the  optimization  problem. 

To  solve  the  optimization  problem  (3.2)-(3.9),  each  optimization  iteration  may 
be  expensive,  because  the  constraints  include  the  full  model.  If  an  output  error 
bound  [60-62,84]  exists,  it  could  be  used  as  the  cost  functional  instead  of  the  true 
output  error.  In  that  case,  the  constraints  only  comprise  the  reduced  model  and  the 
bound  constraints.  As  a  result,  solving  the  optimization  problem  in  this  case  is  much 
less  expensive  since  it  involves  no  full  system  solves.  However,  for  a  general  problem, 
an  error  bound  may  not  be  available.  Alternatively,  an  error  indicator,  for  example 
the  square  of  the  2-norm  of  the  residual,  ||R(<E*xr,  <hxr,  z,  u,  t) |||,  could  be  employed 
(note  that  for  problems  that  result  from  spatial  discretization  of  a  set  of  PDEs  one 
can  use  the  dual  norm  as  in  Section  4.3.1  or  any  weighted-residual  forms).  In  such 
cases,  denote  the  output  error  bound  or  the  norm  of  the  residual  as  <2(x.r,  z,  u,  f);  the 
optimization  problem  now  reads 


max  Q  —  Q(xr,  z,  u,  t)  (3.10) 

x,xr,z 
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subject  to 


At1  R(<E*xr,  <E»xr,  z,  u(f),  t)  = 
xr(0)  = 

^min  —  Z 


0, 

(3.11) 

\hTx°, 

(3.12) 

7 

_  " max  • 

(3.13) 

From  now  on  to  the  end  of  this  chapter,  theoretical  results  will  be  developed  based 
on  the  optimization  problem  using  the  true  output  error,  (3.2)-(3.9).  However,  by 
removing  the  full  model  constraint  and  using  Q(xr,z,u,  t)  instead  of  illy-yrllo  in 
the  cost  functional,  the  results  also  hold  for  the  optimization  problem  (3.10)-(3.13). 

We  denote  the  parameter  vector  that  solves  the  maximization  problem  (3.2)-(3.9) 
by  z*.  Next,  we  compute  the  solution  x(z *,t)  of  the  full  system  at  the  worst-case 
parameter  value  z*.  This  solution  information  is  added  to  the  basis  for  example 
using  the  POD  (note  that  once  the  sample  location  has  been  found,  other  model 
reduction  methods  could  also  be  employed).  The  procedure  is  then  repeated  by 
solving  the  optimization  problem  (3.2)-(3.9)  with  the  updated  basis  <f>.  Thus,  we  are 
using  a  systematic,  adaptive  error  metric  based  on  the  ability  of  the  reduced-order 
model  to  capture  the  outputs  of  interest  in  order  to  choose  the  snapshot  locations  that 
are  locally  the  worst  case  scenarios.  This  adaptive  sampling  approach  is  summarized 
in  the  following  algorithm. 

Algorithm  3.1  Model-Constrained  Adaptive  Sampling  Procedure 

1.  Given  a  reduced  basis  $  and  initial  guess  z°,  solve  the  optimization  problem 
(3.2) -(3.9)  to  find  the  location  in  parameter  space  at  which  the  error  is  maxi¬ 
mized,  i.e.  fiiid  z*  =  arg  max  Q(z). 

2.  If  G( z*)  <  e,  where  e  is  the  desired  level  of  accuracy,  then  terminate  the  algo¬ 
rithm.  If  not,  go  to  the  next  step. 

3.  With  z  =  z* ,  solve  the  full  system  (3.3)  to  compute  the  state  solutions  x(z *,t),  t  = 
(0 ,tf), where  tf  is  some  time  horizon  of  interest.  Use  the  span  of  these  state  so¬ 
lutions  to  update  the  basis  $  such  that  G(z*)  <  e.  Go  to  Step  1. 
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3.1.3  Solution  of  the  Greedy  Optimization  Problem 

In  order  to  reduce  the  offline  cost  of  deriving  the  reduced  model,  we  would  like  to 
minimize  the  cost  of  solving  the  greedy  optimization  problem,  especially  when  large- 
scale  system  equations  appear  as  constraints  in  (3.3).  Therefore,  it  is  important  to 
use  an  efficient  optimization  algorithm  that  allows  us  to  exploit  the  structure  of  the 
system.  In  order  to  solve  the  constrained  optimization  problem  (3.2) — (3.9) ,  we  choose 
to  solve  an  equivalent  bound-constrained  optimization  problem  in  the  z  variables  by 
eliminating  the  state  variables  x  and  xr.  That  is,  we  replace  maxXXrZ  Q (x,  x,r,  z)  with 
maxz  C/(x(z),  xr(z),  z)  =  maxz  Q(z),  where  the  dependence  of  x  and  xr  on  z  is  implicit 
through  the  full  equation  (3.3)  and  reduced  state  equation  (3.6)  (of  course  we  assume 
that  the  full  and  the  reduced  equations  are  well  defined  in  the  sense  that  given  a 
parameter  vector  z  we  can  solve  for  the  full  and  the  reduced  states).  Explicitly,  the 
bound  constrained  optimization  reads 

maxt/(z)  (3-14) 

Z 

subject  to 


^ min  —  ^  —  Z max • 


(3.15) 


While  the  bound-constrained  optimization  problem  (3.14)-(3.15)  is  small,  in  the  sense 
that  the  number  of  optimization  variables  is  small  relative  to  the  size  of  the  state, 
it  may  still  be  expensive  to  solve.  This  is  because  solution  of  the  large-scale  system 
(3.3)-(3.5)  is  still  required  at  each  iteration  of  the  optimization  solver. 

In  the  literature,  there  are  several  methods  to  solve  the  above  bound-constrained 
optimization  problem  [92-101].  In  particular,  the  method  of  Coleman-Li  [98,99],  and 
its  extensions  [95-97, 100, 101]  are  used  here.  Since  our  main  goal  is  to  make  the 
cost  of  solving  the  optimization  problem  as  small  as  possible,  in  the  following  we  will 
combine  the  modified  Coleman-Li  scaling  developed  by  Heinkenschloss  et  al.  [95]  and 
the  subspace  trust  region  interior  reflective  method  in  [101]. 
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The  first  order  necessary  optimality  system  for  the  bound  constrained  problem 
(3.14)-(3.15)  can  be  written  in  the  following  form 

D(z)V£(z)  =  0,  (3.16) 

where  the  diagonal  elements  of  the  Coleman-Li  diagonal  scaling  matrix  D  are  given 
by 

Zi  -  zlin  if  >  0 

Du( z)  =  D?L(z)  =  <  zlax  -  *  if  VQi  <  0  ,  (3.17) 

k  niin{bj  -  zlmml  zlmax  -  zt  j  otherwise 

and  VQi  denotes  the  gradient  of  Q  with  respect  to  z{. 

Next  the  Newton  step  for  (3.16)  at  a  current  optimization  point,  zk,  can  be  written 
as 

M(zfc)s  =  -D(zfc)V£(zfe),  (3.18) 

where  s  is  the  Newton  step,  and  M  =  DV2t?  +  diag(|V(/j|).  Heinkenschloss  et  al.  [95] 
show  that  the  Coleman-Li  scaling  yields  linear  convergence  for  degenerate  cases,  i.e. 
\VGi\  =  0  if  z*  G  {zonin’  zmaX}-  T°  overcome  this,  they  propose  modified  Coleman-Li 
scalings  given  as 

r 

DfiL{z)  if  |Vft|  <  min{zj  -  zlmin,  z. lmax  -  Zi}p  or 
D«(z)  =  *  if  min Ui  -  z*min,  zinax  -  Zi}  <  \WQi\p  (3.19) 

1  otherwise 

for  some  p  >  1. 

The  Coleman-Li  scaling  approach  enables  us  to  use  the  subspace  trust  region 
interior  reflective  Newton  framework,  proposed  in  Branch  et  al.  [101],  to  solve  the 
resulting  bound-constrained  optimization  problem  efficiently.  The  subspace  trust 
region  subproblem  we  need  to  solve  is  given  as 

min  {</9fc(s)  :  ||Ds||2  <  Ak,  s  G  Sk},  (3.20) 

se]Rd 
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where  d  is  the  number  of  parameters,  A*,  is  the  current  trust  region  radius  whose 
updating  rule  can  be  found  in  [101],  and  the  merit  function  </?*,( s)  is  given  as 

(Pk(s)  =  stVQ  +  -s1  D_1Ms.  (3-21) 

We  use  the  conjugate  gradient  (CG)  method  to  determine  the  subspace  S}:  in  which 
the  inexact  Newton  step  s  is  found.  We  terminate  the  CG  subroutine  when  either 
of  the  two  following  conditions  is  satisfied:  (1)  a  negative  curvature  direction  is  en¬ 
countered;  or  (2)  the  norm  of  the  residual  of  the  Newton  system  is  brought  down 
to  a  sufficiently  small  value  relative  to  the  norm  of  the  gradient.  The  subspace  S k 
is  then  constructed  from  the  gradient  and  the  output  of  the  CG  solver.  Finally,  the 
subspace  trust  region  subproblem  (3.20)  is  solved.  This  method  combines  the  rapid 
locally- quadratic  convergence  rate  properties  of  Newton’s  method,  the  effectiveness 
of  trust  region  globalization  for  treating  ill-conditioned  problems,  and  the  Eisenstat- 
Walker  idea  of  preventing  oversolving  [102].  In  addition,  this  method  has  been  nu¬ 
merically  demonstrated  [101]  to  lead  to  faster  convergence  since  it  better  captures 
the  negative  curvature  information  than  does  the  conventional  inexact  Newton- CG 
method  [103,104], 

The  gradient  of  Q  with  respect  to  z,  as  required  by  the  trust  region  subproblem, 
can  be  computed  efficiently  by  an  adjoint  method  which  will  be  developed  in  detail 
for  the  steady  and  unsteady  problems  in  the  following  subsections.  The  Hessian- 
vector  product  as  required  by  CG  is  computed  on-the-fly;  because  it  is  a  directional 
derivative  of  the  gradient  its  computation  similarly  involves  solution  of  state-like  and 
adjoint-like  equations.  Therefore,  the  optimization  algorithm  requires  solution  of  a 
pair  of  state  and  adjoint  systems  at  each  CG  iteration.  Below  is  the  summary  of  the 
optimization  algorithm. 

Algorithm  3.2  Bound- Constrained  Optimization  Solver 

1.  At  the  current  Newton  step  zk ,  compute  the  gradient  V£7( zk)  (requires  a  pair  of 
full  forward  and  full  adjoint  solves,  and  one  pair  of  reduced  forward  and  reduced 
adjoint  solves). 
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2.  Compute  the  subspace  Sk  (each  CG  iteration  requires  a  pair  of  full  forward-  and 
adjoint-like  solves,  and  a  pair  of  reduced  forward-  and  adjoint-like  solves). 

3.  Solve  the  trust  region  subproblem  (3.20). 


4-  Compute  the  following  quantity 


G{  zk  +  s)  -  g(zk) 
<Pk(a) 


(3.22) 


5.  If  p  >  v,  where  v  <  1,  then  set  zk+1  =  zk  +  s  and  A^+i  =  2A^.  Otherwise,  set 
zk+1  =  zk  and  A^+i  =  0.5Afc.  Go  to  Step  2. 


Note  that  a  very  simple  updating  rule  for  zk  and  A^  is  given  in  Step  5  of  Algorithm  3.2 
for  simplicity  of  the  exposition.  In  practice,  a  sophisticated  updating  rule  similar  to 
that  of  Branch  et  al.  [101]  is  used.  Algorithm  3.2  is  repeated  until  some  stopping 
criterion  is  met,  for  example,  when  the  norm  of  the  scaled  gradient  || D || 2  is  less 
than  some  given  tolerance  r. 

Since  the  system  dependence  on  the  parameter  z  is  nonlinear,  in  the  general  case 
the  optimization  problem  (3.2)-(3.9)  is  non-convex.  In  particular,  as  the  greedy 
algorithm  progresses  we  expect  the  cost  functional  to  become  increasingly  multi¬ 
modal,  since  the  error  function  will  be  close  to  zero  (below  the  tolerance  e)  at  each 
of  the  previous  parameter  sample  locations.  It  should  be  noted  that,  while  finding 
the  global  maximum  is  obviously  preferred,  convergence  to  a  local  maximum  is  not 
necessarily  an  adverse  result.  Solving  the  greedy  optimization  problem  is  a  heuristic 
to  systematically  find  “good”  sample  points;  at  a  local  maximum  the  error  is  (locally) 
large.  The  stopping  criterion  applied  in  Step  2  of  Algorithm  3.1  monitors  G(z*),  the 
reduced  model  error  associated  with  the  optimal  solution  z*.  It  is  important  to  note 
that  if  Q( z*)  falls  below  the  desired  error  level,  this  guarantees  only  that  the  local 
error  between  full  and  reduced  model  is  sufficiently  small.  Due  to  the  non-convexity 
of  the  optimization  problem,  it  is  possible  that  larger  errors  may  exist  elsewhere  in 
the  parameter  space. 


57 


3.1.4  Computational  Complexity 


In  this  section,  we  will  approximate  the  cost,  in  terms  of  flop  counts,  of  the  model- 
constrained  adaptive  sampling  procedure,  Algorithm  3.1.  In  order  to  compute  the 
cost,  it  is  helpful  to  point  out  that  the  adaptive  sampling  procedure  involves  two 
main  loops.  The  outer  loop  runs  through  the  number  of  greedy  cycles,  for  which 
Algorithm  3.1  shows  one  cycle.  The  inner  loop  runs  through  the  number  of  Newton 
steps,  one  instance  of  which  is  given  by  Algorithm  3.2.  The  CG  loop,  which  is  inside 
the  inner  loop  occurring  in  Step  2  of  Algorithm  3.2,  is  an  additional  loop  which  is  not 
shown  here.  The  cost  estimation  will  be  only  given  for  steady  problems  of  the  form 
(2.10),  but  the  extensions  to  problems  that  are  nonlinear  in  the  state  vector,  and  to 
unsteady  problems  are  straightforward.  We  further  assume  that  only  the  matrix  A 
depends  on  the  parameters  z  and  is  given  by  the  following  affine  decomposition 

n@ 

A(z)  =  ]Tai(z)At,  (3.23) 

2—1 

where  A,  does  not  depend  on  z  and  0,  is  some  scalar  function  of  z.  Finally,  all  the 
trial  iterates,  zk  +  s,  are  assumed  to  satisfy  p  >  v  in  Step  5  of  Algorithm  3.2.  (In 
the  case  that  this  condition  is  not  satisfied,  the  cost  will  be  slightly  increased  due  to 
additional  Hessian- vector  products  in  Step  2  and  additional  solves  in  Step  4,  because 
one  has  to  restart  Step  2  of  Algorithm  3.2  with  smaller  trust  region  radius.) 

The  order  of  approximating  the  cost  is  as  follows.  The  cost  of  each  Newton  step 
is  first  approximated,  which  is  then  summed  over  all  Newton  steps  to  form  the  cost 
for  each  greedy  cycle.  The  total  offline  cost  is  then  the  sum  of  the  cost  of  all  greedy 
cycles. 

In  order  to  solve  the  full  model,  we  use  the  state-of-the-art  direct  sparse  solver 
UMFPACK  [105-108].  Since  we  consider  problems  that  are  linear  in  the  state  vector, 
we  can  perform  the  LU  factorization  of  the  full  matrix  at  the  beginning  of  Step 
1  of  Algorithm  3.2,  and  then  use  the  full  LU  factors  for  the  full  forward  and  full 
adjoint  solves  and  all  subsequent  full  forward-  and  adjoint-like  solves  required  by  the 
CG  method  within  that  Newton  step.  Similarly,  the  LLI  factorization  for  the  reduced 


model  needs  to  be  computed  once  at  the  beginning  of  Step  1  of  Algorithm  3.2,  and  the 
LU  factors  can  be  then  used  for  the  reduced  forward  and  reduced  adjoint  solves  and 
all  subsequent  reduced  forward-  and  adjoint-like  solves  required  by  the  CG  method. 

Table  3.1  shows  the  cost  in  terms  of  flop  counts  for  each  instance  of  Algorithm  3.2 
(each  Newton  step).  At  the  beginning  of  Step  1  of  Algorithm  3.2,  we  need  one  full  LU 
factorization,  and  its  cost  is  O(LU^),  where  LUj  is  the  number  of  multiplicative  flops 
(*  and  /)  [109].  The  number  of  triangular  solves  for  each  forward  (or  each  adjoint 
solve)  is  one  (it  is  rq,  the  number  of  timesteps,  for  unsteady  problems).  At  the  end 
of  Step  1,  we  need  two  (one  forward  and  one  adjoint)  triangular  solves.  Denote  the 
number  of  Hessian-vector  products  in  Step  2  as  tihv  (which  is  at  most  d  for  the  CG 
method  in  exact  arithmetic)  which  is  assumed  to  be  the  same  for  all  instances  of 
Algorithm  3.2.  Since  each  CG  iteration  requires  one  pair  of  full  forward-  and  adjoint- 
like  solves,  the  number  of  triangular  solves  for  the  CG  solver  is  2 uhv  As  a  result, 
we  have  a  total  of  2 (uhv  +  1)  triangular  solves  for  one  instance  of  Algorithm  3.2;  the 
total  cost  for  all  triangular  solves  is  therefore  2  (nnv  +  1  )0(LU[),  where  LU[  is  the 
number  of  flop  counts  for  each  triangular  solve.  Similarly,  the  flop  counts  related  to 
the  reduced  model  are  shown  in  Table  3.1.  Since  the  reduced  model  is  dense,  the 
explicit  flop  counts  for  LLI  factorization  and  triangular  solve  are  available  in  terms  of 
the  reduced  model  size  rrij,  where  j  denotes  the  jth  greedy  cycle. 

Table  3.1:  Cost  in  terms  of  flop  counts  for  the  bound-constrained  optimization  solver 
for  each  Newton  step  for  steady  problems.  The  true  error  is  used  as  the  cost  functional. 


LLI-factorization 

Triangular  solves 

Number 

Cost 

Number 

Cost 

Full 

1 

0{LUbf) 

2  (riHv  +  1) 

ow) 

Reduced 

1 

O^rrif) 

2  {tlhv  +  1) 

0(mj) 

The  numerical  experiments  in  Sections  4.3.2  and  6.2  suggest  that  the  number  of 
Newton  steps,  and  hence  the  number  of  instances  of  Algorithm  3.2,  scales  linearly 
with  the  number  of  parameters;  we  therefore  assume  the  number  of  Newton  steps  to 
be  O(d).  As  a  result,  the  offline  cost  (flop  counts)  of  computing  a  reduced  basis  with 
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m  basis  vectors  after  nG  greedy  cycles,  excluding  the  cost  of  forming  the  reduced 
model  and  the  cost  of  computing  the  reduced  basis  at  each  greedy  cycle,  scales  as 


/ 


Cost 


TE 


o 


\ 


nG 


d  f  — Wj  “1“  2(1  +  nHv)mj  +  duc^LU ^  +  2(1  +  tihv)LUs  ] 


3= 1 


V 


cost* 


(3.24) 


cost* 


/ 


where  cost*  and  cost*  denote  the  cost  incurred  by  solving  the  intermediate  reduced 
model  and  by  solving  the  full  model,  respectively,  during  the  adaptive  sampling  pro¬ 
cess. 

The  offline  cost  in  terms  of  flop  counts  in  (3.24)  is  for  the  model-constrained 
adaptive  sampling  method  with  the  greedy  optimization  problem  of  the  form  (3.2)- 
(3.9),  that  is,  the  true  error  is  used  as  the  sampling  selection  criteria.  For  the  adaptive 
sampling  method  with  an  error  indicator,  such  as  the  residual  norm,  the  offline  cost 
is  given  by 


Cost 


IE 


O 


\ 


^42 


d  V  p3  +  2(1  +  nHv)j2  +  nG(LU f  +  LUf) 

j= 1  ^  ^  ^ 

.  ^  v  ✓  costF 

\  COStF  / 


(3.25) 


which  involves  only  nG  full  system  solves  at  nG  worst-case  sampling  points.  Since 
only  one  snapshot  at  the  optimal  point  is  computed  in  this  case,  only  one  basis  vector 
is  added  to  the  current  reduced  basis.  As  a  result,  rn3  is  equal  to  j  at  the  end  of  the 
jth  greedy  cycle.  It  can  be  seen  that  if  the  costs  of  the  full  LU  factorization  and  full 
triangular  solves  dominate  the  other  costs,  the  error-indicator  approach  is  much  less 
expensive  than  the  true-error  approach  per  greedy  cycle. 

If  a  good  preconditioner  for  the  CG  solver  is  available  so  that  the  number  of 
Hessian-vector  products  nnv  does  not  depend  on  the  number  of  parameters  d,  the 
offline  cost  of  constructing  reduced  basis  vectors  in  each  greedy  cycle  scales  linearly 
with  the  dimension  of  the  parameter  space,  d.  This  is  clearly  important  as  the  dimen¬ 
sion  of  the  parameter  space  increases.  If,  in  addition,  the  size  of  the  reduced  basis  is 
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constrained  to  be  fixed,  and  hence  the  number  of  greedy  cycles,  tig ,  is  fixed,  it  can  be 
also  proved  that  the  total  offline  cost  of  constructing  the  reduced  basis  scales  linearly 
with  the  dimension  of  the  parameter  space.  However,  if  uq  is  allowed  to  increase 
until  Algorithm  3.1  terminates,  it  will  likely  depend  on  the  dimension  of  the  param¬ 
eter  space,  i.e.  it  might  be  expected  that  larger  input  spaces  require  more  sample 
points,  and  therefore  larger  reduced  basis  size,  as  well  on  complexity  of  dependence  of 
outputs  on  parameters.  In  that  case,  the  total  offline  cost  of  constructing  the  reduced 
basis  no  longer  scales  linearly  with  the  dimension  of  the  parameter  space. 

It  should  be  pointed  out  that  in  general  there  exists  no  good  explicit  approx¬ 
imation  for  0{LUj)  (or  0(LUf))  (of  course  a  very  conservative  upper  bound  is 
0(|n3)  for  0(LUj  )  and  0(2n 2)  for  O(LU^)  from  dense  matrix  theory).  How¬ 
ever,  for  any  well-shaped  finite  element  mesh  [109,  110],  the  flop  counts  can  be 
approximated  as  O(LU^)  +  0(LU[)  =  0(n 3//2)  for  two-dimensional  problems  and 
0(LUj)  +  0{LUg)  =  0(n2)  for  three-dimensional  problems.  Furthermore,  the  stor¬ 
age  requirement  for  the  LU  factors  is  0(nlogn )  for  two-dimensional  problems  and 
(9(n4/3)  for  three-dimensional  problems.  Therefore,  for  three-dimensional  problems,  a 
direct  solver  may  not  be  possible,  and  an  iterative  sparse  solver  is  an  alternative  [111]. 

To  estimate  the  cost  of  forming  the  reduced  model  after  each  adaptive  cycle,  we 
first  consider  the  error-indicator  approach  for  simplicity.  The  discussion  for  the  cost 
estimation  of  the  true-error  will  then  follow.  Since  the  matrix  A  is  assumed  to  be  of 
the  form  (3.23),  we  compute  the  cost  of  forming  the  reduced  model  explicitly.  From 
the  reduced  equation  (2.15),  the  reduced  matrix  is  given  as 

ne  ne 

Ar  =  0fc(z)0i(z)$rA^A^,  (3.26) 

k= 1  /=1 

which  suggests  that  we  need  to  compute  and  store  the  matrices 

Af;1  =  $tA[A;$.  (3.27) 

It  should  be  noted  that  the  dimension  of  the  reduced  matrix  Af  computed  at  the 
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end  of  the  (j  +  l)th  cycle  is  (j  +  1)  x  (j  +  1),  but  its  j  x  j  principle  submatrix  is 
the  reduced  matrix  at  the  end  of  the  jth  cycle.  Therefore,  one  only  needs  to  update 
the  (j  +  l)th  row  and  (j  +  l)th  column  of  after  the  (j  +  l)th  greedy  cycle.  Next, 
denote  nnz  to  be  the  number  of  nonzeros  of  A©  and  assume  it  is  the  same  for  all  A*. 
To  compute  the  first  j  elements  of  the  ( j  +  l)th  column,  we  need  to  compute 

$TA^$i+1,  (3.28) 

where  <3>,)+i  is  the  new  basis  vector  computed  at  the  end  of  the  (j  +  l)th  cycle.  Multi¬ 
plying  from  right  to  left,  the  cost  is  approximately  given  as  2jn  +  2nnz,  which  is  also 
the  cost  of  computing  the  first  j  elements  of  the  (j  +  l)th  row.  Finally,  the  cost  of 
computing  the  (j  +  l)th  element  of  the  (j  +  l)th  column  is  given  as  2 n  +  2 nnz.  The 
total  cost  of  updating  the  reduced  matrix  A^}  is  given  as  2(2 j  +  l)n  +  6nnz.  Therefore 
the  total  cost  to  compute  all  the  reduced  matrices  at  the  end  of  the  (j  +  l)th  greedy 
cycle  is  given  as 

n@[2(2j  +  1  )n  +  6  nnz\,  (3.29) 

and  the  corresponding  storage  requirement  is  (j  + 1)2.  It  can  be  seen  that,  depend¬ 
ing  on  n© ,  j,  n,  and  nnz,  the  cost  of  computing  the  intermediate  reduced  matrices 
can  form  a  considerable  portion  of  the  total  offline  cost. 

For  the  true-error  case,  the  cost  of  forming  the  reduced  model  after  each  greedy 
cycle  is  the  same  as  that  of  the  error-indicator  case  if  only  the  snapshot  at  the  optimal 
parameter  point  is  used  to  update  the  reduced  basis.  However,  as  will  be  discussed  in 
Section  3.1.6,  the  size  of  the  reduced  model  depends  on  how  we  update  the  reduced 
basis,  and  hence  the  cost  could  be  slightly  increased.  Nonetheless,  the  procedure  for 
estimating  the  cost  of  forming  the  reduced  basis  is  similar. 

It  must  be  also  mentioned  that  computing  reduced  basis  vectors  and  using  them 
to  update  the  reduced  basis  can  take  a  considerable  amount  of  time  (due  to  dense 
linear  algebra  operations  on  these  dimension  n  reduced  basis  vectors)  and  storage, 
depending  on  how  many  basis  vectors  are  generated  and  on  the  size  of  the  full  model, 
n. 
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3.1.5  Initialization 


Initially,  there  are  no  basis  vectors  in  the  reduced  basis;  it  is  therefore  natural  to 
choose  the  initial  basis  as  the  empty  set,  =  0,  and  the  reduced  model  is  a  zero- 
order  approximation  of  the  full  model. 

It  should  be  also  pointed  out  that  the  Algorithm  3.2  starts  with  an  initial  guess 
z°  in  the  parameter  space  and  moves  iteratively  towards  a  local  maximizer.  To 
avoid  convergence  to  a  local  maximum  close  to  a  previous  sample  location,  and  hence 
to  explore  the  parameter  space  better,  a  random  initialization  of  the  optimization 
variables  z  is  used.  An  initial  guess  is  accepted  only  if  it  is  “far  enough”  away  from 
the  previous  sample  locations  and  its  corresponding  cost  functional  is  larger  than 
e.  In  particular,  the  smallest  allowable  distance  between  an  initial  guess  and  all  the 
existing  sample  locations  is  chosen,  for  example,  to  be  0.5  min i{zlmax  —  zlminj. 

3.1.6  Updating  the  Reduced  Basis 

Recall  that  the  purpose  of  Algorithm  3.1  is  to  sample  the  (locally)  optimal  parameter 
points  in  the  parameter  space.  Therefore,  the  snapshots  at  these  optimal  sample 
locations  arc  important,  and  it  is  natural  to  use  these  snapshots  to  update  the  reduced 
basis  as  proposed  in  Step  3  of  Algorithm  3.1.  However,  if  the  true  output  error  is  used 
as  the  cost  functional,  besides  the  snapshots  computed  at  the  optimal  solutions,  the 
solution  snapshots  and  the  adjoint  solutions  at  each  Newton  step  are  also  available 
as  part  of  the  optimization  process;  although  these  snapshots  may  not  add  much 
information  as  those  at  the  optimal  parameters,  they  can  be  used  to  improve  the 
reduced  basis.  Given  all  these  snapshots,  some  approaches  to  update  the  reduced 
basis  are  given  as  follows 

1.  Use  only  the  snapshots  at  the  optimal  parameter  points,  z*,  to  update  the 
reduced  basis,  for  example,  using  the  Gram-Schmidt  procedure. 

2.  Store  all  snapshots,  and  then  perform  the  POD  method  on  the  complete  snap¬ 
shot  set,  which  comprises  those  of  the  current  greedy  cycle  and  all  previous  ones. 
For  large-scale  problems,  this  approach  is,  however,  expensive  both  in  terms  of 
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storage  and  in  terms  of  computing  time.  A  less  expensive  approach  is  first  to 
perform  the  POD  method  on  just  the  snapshots  of  the  current  greedy  cycle  in 
order  to  extract  dominant  POD  basis  vectors.  These  dominant  POD  vectors  are 
then  added  to  the  current  reduced  basis  using  Gram-Schmidt  orthogonalization. 
Another  approach  is  to  perform  a  second  POD  computation  on  the  snapshot  set 
comprising  the  current  reduced  basis  and  the  newly  computed  POD  vectors.  In 
that  case,  since  the  POD  basis  vectors  are  of  unit  length,  an  appropriate  scaling, 
for  example  using  the  corresponding  singular  values,  needs  to  be  done  before  the 
second  POD  computation.  Otherwise,  the  important  information  compressed 
in  the  dominant  POD  vectors  of  the  current  reduced  basis  and  of  the  POD  basis 
of  the  current  greedy  cycle,  may  be  lost,  resulting  in  a  reduced  basis  that  is  poor 
in  quality. 


3.  Add  all  snapshots  of  the  current  greedy  cycle  to  the  current  reduced  basis 
using  Gram-Schmidt  orthogonalization.  However,  care  must  be  taken  since  this 
method  potentially  makes  both  the  offline  and  the  online  stages  expensive  due 
to  a  large  number  of  vectors  in  the  reduced  basis.  Yet,  the  reduced  model  may 
not  be  accurate  because  not  all  the  snapshots  are  important. 


4.  Solve  an  (inner)  optimization  problem  to  find  the  basis  that  minimizes  the 
output  error  at  the  sample  points  at  which  the  snapshots  are  computed  [112]. 


If  the  error  indicator  is  used,  only  snapshots  at  the  optimal  parameter  points 
are  computed.  As  a  result,  if  the  number  of  optimal  parameter  points,  and  hence 
the  number  of  snapshots,  is  small,  one  can  simply  use  the  first  updating  approach, 
as  discussed  above.  However,  if  the  number  of  snapshots  is  large,  all  the  above 
approaches  can  be  employed. 
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3.2  An  Analysis  of  the  Adaptive  Sampling  Ap¬ 


proach 

We  have  proposed  an  adaptive  sampling  method  for  parameter-dependent  problems 
in  Section  3.1.  Recall  that  the  optimization  problem  in  each  adaptive  sampling  pro¬ 
cedure  is  a  PDE-constrained  optimization  problem.  As  a  result,  the  question  arises 
under  what  conditions  the  adaptive  sampling  procedure  works.  In  particular,  we  need 
to  answer  the  following  questions:  1)  Does  the  PDE-constrained  optimization  prob¬ 
lem  in  each  adaptive  cycle  have  a  solution?  2)  If  it  does,  is  the  solution  unique?  That 
is,  the  existence  and  uniqueness  of  the  solution  of  the  optimization  problem  need  to 
be  addressed. 

3.2.1  Existence  and  Uniqueness  of  a  Maximizer  in  Each  Adap¬ 
tive  Cycle 

In  this  section,  we  prove  that  a  solution  of  the  optimization  problem  (3.14) — (3.15),  and 
hence  the  optimization  problem  (3.2)-(3.9)  or  (3.10)— (3.13),  exists,  and  discuss  the 
uniqueness  of  that  solution.  To  begin,  let  us  recall  one  of  the  fundamental  theorems 
about  continuous  functions. 

Theorem  3.1  If  Q{z\, . . . ,  zf)  :  Q  C  IRd  — >  IR  is  continuous  and  Vt  is  a  compact 
subset  of  Md ,  then  there  is  a  maximum  point  z  efl  such  that  Q(z)  <  Q(  z),Vz  G  D. 

Proof:  The  proof  can  be  found,  for  example,  in  [113].  □ 

This  theorem  shows  that  the  greedy  optimization  problem  in  each  adaptive  cycle 
has  at  least  one  solution  by  the  following  corollary. 

Corollary  3.1  In  each  adaptive  cycle,  assume  that  the  cost  functional  is  a  continuous 
function  in  the  parameter  z,  then  there  exists  at  least  one  solution  for  the  optimization 
problem  (3.14)-(3.15). 

Proof:  Denote  fl  to  be  the  parameter  set  defined  by  the  bound  constraints  (this  is 
called  a  d-cell  in  mathematics  [113]).  Then,  it  can  be  proved  that  a  d-cell  is  compact 
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(the  proof  requires  a  few  other  theorems  in  set  topology  study  and  therefore  we  simply 
use  this  result). 

On  the  other  hand,  the  states  are  eliminated  so  that  the  cost  function  G(z)  is 
only  a  function  of  the  parameter  z.  Therefore,  if  cost  functional  Q( z)  is  a  continuous 
function  on  0,  Theorem  3.1  applies  and  Corollary  3.1  is  proved.  That  is,  there  exists 
a  solution  (a  global  maximizer  according  to  the  theorem)  to  the  optimization  problem 
(3.14)-(3.15)  in  each  adaptive  cycle.  □ 

Clearly  uniqueness  is  not  guaranteed  in  the  general  case  since  there  could  be  many 
global  maximizers. 

3.2.2  Properties  of  the  Adaptive  Sampling  Approach 

Next,  some  important  properties  of  the  adaptive  sampling  approach  will  be  discussed. 

Theorem  3.2  Assume  that  the  full  model  is  linear  in  state  x.  Then,  in  the  kth 
adaptive  cycle,  the  cost  functional  is  less  than  e  at  all  the  maximizers  found  in  the 
previous  cycles  k  <  k. 

Proof:  Recall  that  in  Step  3  of  Algorithm  3.1  the  span  of  the  state  solutions  at  the 
local  maximizers  found  in  previous  cycles  k  <  k  are  used  as  basis  vectors  such  that 
the  cost  functional  at  these  local  maximizers  is  less  than  e.  Furthermore,  we  proved 
in  Chapter  2  that,  for  problems  that  are  linear  in  state  x,  as  the  reduced  basis  is 
enriched,  the  reduced  model  error  cannot  increase.  As  a  result,  the  cost  functional  in 
the  kth  adaptive  cycle  is  less  than  e  at  all  the  maximizers  found  in  the  previous  cycles 
k  <  k.  □ 

As  a  consequence  of  the  above  theorem,  the  below  corollary  is  an  important  result 
for  the  adaptive  sampling  approach. 

Corollary  3.2  Assume  that  the  full  model  is  linear  in  state  x.  Then,  the  adaptive 
sampling  approach  will  never  sample  at  the  previous  sampled  points  in  the  parameter 
space. 
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Proof:  By  definition  in  (3.2),  the  cost  functional  is  non- negative.  To  prove  the 
corollary,  it  is  sufficient  to  show  that  in  the  kth  adaptive  cycle  the  maximizer  must  be 
different  from  the  maximizers  found  in  the  previous  cycles.  First,  recall  that  the  cost 
functional  in  the  current  adaptive  cycle  is  smaller  than  e  at  all  previous  maximizers, 
as  proved  in  Theorem  3.2.  Second,  we  only  start  at  an  initial  guess  where  the  cost 
function  is  greater  than  e.  Third,  the  optimization  solver  only  accepts  an  iterate 
if  the  cost  functional  is  larger  than  that  at  the  previous  iterate.  Using  these  three 
facts  we  conclude  that  the  cost  functional  at  a  new  maximizer  must  be  larger  than  e. 
Therefore,  the  maximizer  found  in  the  kth  adaptive  cycle  must  be  different  from  the 
previous  maximizers.  □ 

Even  though  we  do  not  yet  have  proofs  that  the  results  in  this  section  are  true 
for  general  problems  that  are  nonlinear  in  both  state  and  parameters,  this  is  not  a 
severe  limitation  if,  for  example,  one  uses  the  empirical  interpolation  method  [61,62, 
68,69,114,115]  to  pre-process  the  nonlinear  terms  by  a  linear  combination  of  empirical 
interpolation  basis  functions.  In  that  case,  the  results  in  this  section  still  apply  for 
the  pre-processed  model.  Of  course,  if  one  can  also  prove  that,  for  a  problem  at  hand 
that  is  nonlinear  in  state  x,  the  reduced  model  is  improved  as  the  reduced  basis  is 
enriched,  both  Theorem  3.2  and  Corollary  3.2  also  hold  for  that  nonlinear  problem. 

We  next  derive  the  optimization  problem  in  each  adaptive  cycle  and  its  optimality 
system  for  a  class  of  steady  and  unsteady  problems  that  arc  linear  in  the  state  vector, 
but  depend  nonlinearly  on  the  parameters. 


3.3  Steady  Problems 

Consider  a  general  large-scale  parameter-dependent  steady  problem  that  is  linear  in 
the  state  vector  x 


A(z)x  =  B(z);  y  =  C(z)x, 


<  z  <  z. 


-imin  _  _ 


(3.30) 
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As  discussed  in  Section  2.1,  a  projection-based  model  order  reduction  technique  yields 
the  reduced  system  of  the  form 


A,.(z)xr  Br(z),  yr  Cr(z)xr,  z min  ^  z  ^  zmax, 


(3.31) 


where  A,r(z)  =  th1  A(z)<fr,  Br(z)  =  \l/TB(z)  and  Cr(z)  =  C(z)<f>.  Note  that  the  full 
and  the  reduced  models  were  given  in  (2.10)  and  (2.12);  they  are  repeated  here  for 
clarity.  Next,  define  the  cost  functional 


1.,  ,2  1 


£(x,x,r,z)  =  ^||y-yr||o  =  ^lly-yrlll, 


(3.32) 


which  describes  the  output  error  between  the  full  and  the  reduced  models.  Here, 
||  •  ||2  is  the  Euclidean  2-norm  (any  weighted  combination  of  the  outputs  could  be 
used  instead).  The  greedy  optimization  problem  in  each  adaptive  cycle  then  reads 


max  Q  =  -  ||C(z)x  —  Cr(z)xr||2, 

X,Xr,Z  2 


(3.33) 


subject  to 


A(z)x  =  B(z), 

(3.34) 

A.r(z)xr  =  Br(z), 

(3.35) 

^min  A  ^  A  Zmax. 

(3.36) 

The  optimality  conditions  for  the  constrained  optimization  problem  (3.33) — (3.36) 
can  be  derived  by  defining  the  Lagrangian  functional 

£(x,  xr,  z,  A,  Ar)  =  C/(x,  xr,  z)  +  At  [A(z)x  —  B(z)] 

+  Aj  [Ar(z)xr  -  Br(z)] ,  (3.37) 

where  A  G  Rn  and  Ar  G  Rm  are  the  full  and  the  reduced  adjoint  variables  that  respec¬ 
tively  enforce  the  full  and  the  reduced  equations.  Note  that  the  bound  constraints  are 


excluded  and  treated  separately  as  in  Section  3.1.3.  The  first-order  Karush-Kuhn- 
Tucker  optimality  system  can  be  derived  by  taking  derivatives  of  the  Lagrangian  with 
respect  to  the  adjoints,  states,  and  parameter  as  follows: 


•  Setting  the  derivative  of  the  Lagrangian  with  respect  to  A  to  zero  yields  the  full 
equations  (3.34). 

•  Setting  the  derivative  of  the  Lagrangian  with  respect  to  \r  to  zero  yields  the 
reduced  equations  (3.35). 

•  Setting  the  derivative  of  the  Lagrangian  with  respect  to  x  to  zero  yields  the  full 
adjoint  equations 


Ar(z)A  =  Cr(z)  [Cr(z)xr  —  C(z)x  . 


(3.38) 


Setting  the  derivative  of  the  Lagrangian  with  respect  to  xr  to  zero  yields  the 
reduced  adjoint  equations 


A^(z)Ar  =  CjT(z)  [C(z)x  -  Cr(z)x.r  . 


(3.39) 


Setting  the  derivative  of  the  Lagrangian  with  respect  to  z%  to  zero  yields  the 
optimality  conditions 


\dC,  . 

dCr,  . 

T 

dA 

dB,  J 

dzt  (Z)X  j 

[C(z)x  —  Cr(z)xr]  +  Ar 

k(z,x" 

+  A 


dK,  ,  dBr 

Z  x,r —  z 


=  0. 


d Zi  v  '  dzi 

The  reduced  gradient  of  the  cost  functional  is  then  given  by 


d_Q_ 

d  Zi 


(3.40) 


\dC,  . 

dCr,  . 

T 

fclA.  . 

dB.  J 

[ay(z)x- 

dz,  (z)xj 

[C(z)x  —  Cr(z)xr]  +  AT 

1 

1  O 
_ 1 

+  A I 


dAr  OBr 

"&r(z)Xr  -  sy(z) 


(3.41) 
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where  the  full  state  x,  reduced  state  x.r,  full  adjoint  variable  A  and  reduced  adjoint 
variable  Ar  respectively  satisfy  the  full  forward  equation  (3.34),  the  reduced  forward 
equation  (3.35),  the  full  adjoint  equation  (3.38),  and  the  reduced  adjoint  equation 
(3.39). 

3.4  Unsteady  Problems 

For  clarity,  we  recall  the  results  in  Section  2.1  that  given  a  general  parametrized  LTI 
dynamical  system 


E(z)x 

=  A(z)x+B(z)u, 

(3.42) 

y 

=  C(z)x, 

(3.43) 

x(0) 

(3.44) 

7 

" mm  — 

^  ^  Z inaxi 

(3.45) 

a  projection-based  model  order  reduction  technique  yields  the  reduced  system  of  the 
form 


Er(z)xr  =  A.r(z)xr  +  Br(z)u,  (3.46) 

yr  =  Cr(z)xr,  (3.47) 

xr(0)  =  ^Tx°,  (3.48) 

Z min  —  ^  —  Z maxi  (3.49) 


where  Er(z)  =  VTE(z )$,  Ar(z)  =  ^TA(z)$,  Br(z)  =  ^TB(z),  Cr(z)  =  C(z)$. 
We  define  the  cost  functional 


£(z)  =  ^lly-yrllo  =  \  I  I|y-yrll2 dt  =  \fo  l|C(z)x-cr(z)xr|||di,  (3.50) 

which  describes  the  error  between  the  full  and  reduced  models  over  the  parameter 
space  z,  integrated  over  some  time  horizon  of  interest  tf.  The  greedy  optimization 
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problem  in  each  adaptive  cycle  then  reads 


1  ftf 

max  Q  =  -  /  ||C(z)x  —  Cr(z)xr||9  dt,  (3.51) 

x,xr.z  2  Jo 

subject  to 


E(z)x 

=  A(z)x  +  B(z)u, 

(3.52) 

x(0) 

(3.53) 

E,.(z)xr 

=  Ar(z)xr  +  B.r(z)u, 

(3.54) 

xr(0) 

=  H>Tx°, 

(3.55) 

Z/nin  L 

2  L;  Zmax- 

(3.56) 

The  optimality  conditions  for  the  constrained  optimization  problem  (3.51)-(3.56) 
can  be  derived  by  defining  the  Lagrangian  functional 

1  ftf 

£(z.  X,  xr,  A,  Ar,  7, 7r)  =  -  /  ||C(z)x  -  C.r(z)xr||7df 

*  Jo 

rtf  rtf 

+  /  A7  [E(z)x  —  A(z)x  —  B(z)u]  dt  +  /  A)T  [Er(z)xr  —  A,r(z)x,r  —  Br(z)u]  dt 

W  [x(0)  -  X°]  +  [xr(0)  -  ^ }  ,  °  (3.57) 

where  A  G  1RM  and  Ar  G  !Rm  are  the  full  and  the  reduced  adjoint  variables  that  respec¬ 
tively  enforce  the  full  and  the  reduced  equations.  Two  additional  adjoint  variables 
7  G  IR"  and  7r  G  ]Rm  are  used  to  enforce  the  initial  conditions  for  the  full  and  the  re¬ 
duced  models,  respectively.  Note  that  the  bound  constraints  are  excluded  and  treated 
separately  as  in  Section  3.1.3.  The  first-order  Karush-Kuhn- Tucker  optimality  system 
can  be  derived  by  taking  variations  of  the  Lagrangian  with  respect  to  the  adjoints, 
states,  and  parameter  as  follows: 

•  Setting  the  first  variation  of  the  Lagrangian  with  respect  to  A  to  zero,  and 
arguing  that  the  variation  of  A  is  arbitrary  in  (0, 7/) ,  yield  the  full  equation 
(3.52). 
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Setting  the  first  variation  of  the  Lagrangian  with  respect  to  Ar  to  zero,  and 
arguing  that  the  variation  of  Ar  is  arbitrary  in  (0,  tf),  yield  the  reduced  equation 
(3.54). 


•  Setting  the  derivative  of  the  Lagrangian  with  respect  to  7  to  zero  yields  the  full 
initial  condition  equation  (3.53). 


•  Setting  the  derivative  of  the  Lagrangian  with  respect  to  to  zero  yields  the 
reduced  initial  condition  equation  (3.55). 


•  Setting  the  first  variation  of  the  Lagrangian  with  respect  to  x  to  zero,  and 
arguing  that  the  variation  of  x  is  arbitrary  in  (0,  tf),  at  t  —  0  and  at  t  =  tf, 
yield  the  full  adjoint  equation,  final  condition  and  definition  of  7 


Et(z)A  +  At(z)A 

a  on 

7 


C7  (z)  [C(z)x  —  Cr(z)x.r] , 

0, 

ET(z)A(0). 


(3.58) 

(3.59) 

(3.60) 


•  Setting  the  first  variation  of  the  Lagrangian  with  respect  to  x.r  to  zero,  and 
arguing  that  the  variation  of  x.r  is  arbitrary  in  (0 ,tf),  at  t  —  0  and  at  t  =  tf, 
yield  the  reduced  adjoint  equation,  final  condition  and  definition  of 


E^(z)Ar  +  A)f(z)Ar 


C^(z)  [Cr(z)x.r  -  C(z)x] , 


(3.61) 


A r(T )  —  0, 

7 r  =  E^(z)Ar(0). 


(3.62) 

(3.63) 


•  Setting  the  derivative  of  the  Lagrangian  with  respect  to  zt  to  zero  yields  the 
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optimality  condition 
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The  reduced  gradients  of  the  cost  functional  with  respect  to  z  are  then  given  by 
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where  the  full  state  x,  reduced  state  xr,  full  adjoint  variable  A  and  reduced  variable 
\r  respectively  satisfy  the  full  forward  equation  (3.52),  the  reduced  forward  equa¬ 
tion  (3.54),  the  full  adjoint  equation  (3.58)— (3.59) ,  and  the  reduced  adjoint  equation 
(3.61)— (3.62).  In  addition,  the  adjoints  7  and  qr  satisfy  (3.60)  and  (3.63),  respectively. 

Note  that  derivation  of  the  optimization  problem  for  discrete  time  systems  and 
their  corresponding  optimality  conditions  is  similar,  and  hence  omitted  here. 

A  question  one  may  immediately  ask  is  what  kind  of  input,  u  =  u(t)  we  should 
use  in  the  optimization  problem.  The  answer  is  given  by  the  following  proposition. 


Proposition  3.1  Assume  the  full  system  (3.52)  is  a  single-input  multi-output  (SIMO) 
dynamical  system.  If  the  output  error  norm  for  the  impulse  response  is  bounded,  i.e., 


lls(z)llo  =  /  ll(Mr)  -hr(r))\\l  dr  < 
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where  hit )  and  hr(t )  are  the  impulse  responses  of  the  full  and  reduced  models,  respec¬ 
tively,  then  the  error  norm  for  an  arbitrary  input  u  with  finite  ?iorm  is  also  bounded, 
i.e., 


l|G(z)|| 


o 


u 


<  e 


where  the  norm  of  the  input  is  given  as  ||-u||f 


Iof  Iiu(t-r)2  dr  dt. 


Proof:  Using  the  definition  of  the  output  norm,  convolution  integral  and  Cauchy- 
Schwarz  inequality  gives 
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□ 

It  should  be  pointed  out  that  it  is  sufficient  to  consider  a  SIMO  system  since  a 
multi-input  multi-output  system  is  equivalent  to  a  superposition  of  SIMO  systems. 

The  above  proposition  is  important  for  the  greedy  optimization  process,  that  is, 
we  need  to  only  consider  the  impulse  input  in  the  adaptive  sampling  procedure.  Once 
the  error  of  the  reduced  model  is  small  for  the  impulse  input,  it  will  be  small  for  an 
arbitrary  input  with  finite  norm.  While  the  result  in  Proposition  3.1  is  general,  in 
practice  the  frequency  content  of  the  input  is  sometimes  known  a  priori.  In  that  case, 
an  appropriate  input  (i.e.  not  the  impulse)  with  the  given  frequency  content  could 
be  used. 

For  problems  that  are  linear  in  state  vector  x,  one  can  transform  the  full  and  the 
reduced  equations  from  the  time  domain  into  the  frequency  domain.  In  the  frequency 
domain,  the  problem  can  be  considered  as  a  steady  system,  and  the  forcing  input 
is  replaced  by  the  frequency  u.  The  trade-off  here  is  that  the  size  of  the  system  is 
double,  though  the  system  has  some  structure,  and  the  frequency  u>  is  an  additional 
parameter  which  is  assumed  to  be  in  the  range  w>min  Umax,  where  [w>min^max\ 
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is  the  frequency  range  of  interest.  However,  conceptually  one  could  solve  a  greedy 
optimization  problem,  by  exploiting  the  special  structure  of  the  problem,  to  select 
frequencies  within  this  range  at  which  snapshots  should  be  computed. 
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Chapter  4 


Steady  Thermal  Fin  Heat 
Conduction  Models 


In  this  chapter,  we  consider  the  thermal  fin  design  problem  adopted  from  [116].  While 
the  thermal  fin  geometry  in  [116]  is  fixed,  it  is  allowed  to  vary  in  this  thesis.  As  a 
result,  our  problem  is  described  by  thirty-four  parameters.  The  detail  of  the  problem 
in  both  physical  domain  and  computational  domain  is  first  described.  The  application 
of  the  model-constrained  greedy-based  sampling  approach  developed  in  Chapter  3  to 
obtain  reduced  models  for  the  thermal  hn  problem  is  then  performed  for  different 
number  of  parameters.  Finally,  the  thermal  hn  optimal  design  problem  with  11 
parameters  using  the  reduced  model  is  presented  and  compared  to  that  obtained 
using  the  full  model. 

4.1  Physical  Domain  Formulation 

The  two-dimensional  thermal  hn  is  shown  in  Figure  4-1.  It  consists  of  the  vertical  post 
and  four  horizontal  subhns.  The  purpose  of  the  thermal  hn  is,  for  example,  to  conduct 
heat  from  some  machine,  which  is  attached  to  the  root,  through  the  large-surface- 
area  subhns  to  the  surrounding  flowing  air.  The  vertical  post  has  eight  sub-regions 
denoted  as  09, . . . ,  f216  with  corresponding  thermal  conductivities  Kg, ,  kw  (e.g. 
these  regions  could  be  made  of  different  materials).  Each  subhn  has  two  different 
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sub-regions,  e.g.  f^i-i  and  Cl2i  for  the  ith  subfin,  where  i  =  1,  ...,4,  and  these 
sub-regions  have  different  thermal  conductivities  denoted  as  Ki, . . . ,  k$.  In  addition, 
the  size  of  all  the  sub-regions  of  the  post  and  the  subfins  could  be  varied  and  they 
are  denoted  as  bi, . . .  ,b77,  as  shown  in  Figure  4-1.  Another  parameter  of  interest  is 
the  Biot  number,  Bi,  which  characterizes  the  convective  heat  to  the  air  at  the  fin 
surfaces.  Therefore,  we  have  in  total  thirty-four  parameters,  which  are  represented 
by  the  vector  of  parametric  inputs  z  as  z  =  {z\, . . . ,  Z34},  where  Zi  =  /c*,  i  =  1, . . . ,  16, 
z17  =  Bi ,  and  zl7+j  =  bj,  j  =  1, . . . ,  17. 
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Figure  4-1:  The  thermal  fin  geometry  in  the  physical  domain. 

The  steady-state  temperature  distribution  within  the  fin,  w(z),  is  governed  by  the 
following  elliptic  PDE 


— KiV2wl  =  0  in  Oj,  i  —  1, . . . ,  16, 


(4.1) 


where  wl  denotes  the  restriction  of  w  to  hi*,  V  =  Jri  +  J^j,  and  V2 


a2  ,  a2 
as?  ^  7^’ 
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where  i  and  j  are  the  unit  vectors  pointing  along  the  x—  and  y— direction.  The 
continuity  of  temperature  and  heat  flux  at  the  conductivity-discontinuity  interfaces 
T™1  =  dCli  fi  dClj  for  two  adjacent  regions  0,:  and  Clj,  where  <9f 2*  and  dVLj  denote  the 
boundary  of  O,  and  Qj  respectively,  are  ensured  by  the  following  interface  condition 


w 


wJ 


-Ki{S/wl  ■  n*)  =  —K,j('Vwj-hj) 


on  T 


int 
ij  > 


(4.2) 


where  n*  and  lV  denote  the  outward  normal  of  O,  and  Clj  on  the  interface  T™1, 
respectively.  In  order  to  model  the  convective  heat  losses  on  the  external  surface  of 
a  region  Qt,  i.e.  rfxt  =  d fh  \  T”1*  and  rfxt  ^  0,  we  use  the  following  Robin  boundary 
condition 

-KitW  •  n*)  =  Bi  wi  on  Tf*  if  Ff*  ^  0,  i  =  1, . . . ,  16.  (4.3) 


Finally,  to  model  the  heat  source  at  the  root,  the  Neumann  flux  boundary  condition 
is  imposed  as 

-k9(Vw9  •  n9)  =  -1  on  froot.  (4.4) 


For  this  particular  problem,  the  output  of  interest  is  the  average  temperature  over 
the  whole  thermal  fin 


Ei=i  Int  w  d ^ 
^2i= i  Ify  1 


(4.5) 


Following  Ref.  [116],  it  can  be  showed  that  the  temperature  distribution  solution 
w  belongs  to  the  Hilbert  space  Hl{ 0),  where  0  =  and  satisfies  the  following 

weak  form 

a(w,v )  =  i(v),  Vv  G  Ff1(f2),  (4.6) 

where  the  bilinear  form  a  is  given  as 

16  „  16 
a(w,  v )  =  ^  /  Vw  ■  Vv  dQ,;  +  Bi 

1=1  i=  1 


wv  dr 


ext 


I  pext 


(4.7) 
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and  the  linear  form  l  as 


£(v)  =  /  V  dTmot  (4.8) 

J  proot 


4.2  Computational  Domain  Formulation 

As  a  common  practice,  it  is  useful  to  transform  the  physical  domain  to  a  computa¬ 
tional  (or  reference)  domain  to  carry  out  the  numerical  calculation.  A  clear  advantage 
is  that  one  only  needs  to  mesh  the  computational  domain  once,  and  the  computa¬ 
tional  mesh  can  be  used  to  solve  the  problem  with  any  set  of  thermal  conductivities, 
lengths  and  Biot  number.  The  computational  domain  is  chosen  as  in  Figure  4-2  in 
which  the  dimensions  of  all  computational  regions  are  also  presented. 
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Figure  4-2:  The  thermal  fin  geometry  in  the  computational  domain. 


Using  a  few  simple  linear  transformation  rules  from  xy— coordinates  to  xy— coordinates, 
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one  can  write  the  weak  form  (4.6)  as 


a(w,v )  =  £(v),  Vn  e  ih1(h2),  Q  = 
where  the  linear  form  i  is  now  given  as 


(4.9) 


flv)  =b„  v  dr°°\ 

J  proot 


(4.10) 


and  the  bilinear  form  a  as 


a(w,  v )  = 


i= i  1 


dw  dv  ^  f  dw  dv 

TT-  TT"  ^0,;  +  02j  /  TT-  7T" 

n,  ox  ox  Jo  dy  dy 


dih  1  +  V"'  @7+32  [ 

J  j=l  Jrtf 


wv  drif, 
(4.11) 

where  T^xt  denotes  the  external  boundary  corresponding  to  bj,  and  0,;  are  given  as 
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The  output  of  interest  is  now  given  as 

=  E2i  A  /n.  w  dttt 

v  E‘h  ft  41  <«v 

where  are  given  by 


(4.12) 
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Next  using  the  finite  element  method  with  triangular  linear  elements,  for  example, 
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the  weak  form  (4.9)  can  be  written  in  matrix  form  as 


A(z)x  =  B(z),  y  =  C(z)x, 


(4.13) 


where  the  vector  x  consists  of  all  nodal  temperature  values.  The  matrices  A(z),  B(z), 
and  C(z)  are  as  follows 


A(z) 

48 

=  £0‘(Z)A- 
i  —  1 

(4.14) 

B(z) 

l —  1 

—  ^34B  q, 

(4.15) 

16 

C(z) 

=  EA-(z)c»- 

i= 1 

(4.16) 

where  Aqi,  Bg  and  Cqi  are  the  appropriate  finite  element  matrices  that  do  not  depend 
on  the  parameters  z,  and  A,:  are  given  as 


A*(z) 


_ A(z) _ 

1  d(t> 


(4.17) 


It  should  be  pointed  out  that  if  one  chooses  the  average  temperature  at  the  root  as 
the  output  of  interest,  the  output  matrix  C  is  then  exactly  which  does  not  depend 
on  the  parameters  z. 


The  matrices  in  (4.14)-(4.16)  admit  an  affine  decomposition  in  which  the  affine 
coefficients,  0,  and  A*,  are  nonlinear  functions  of  the  parameters  z.  This  affine  de¬ 
composition  together  with  the  projection-based  model  reduction  technique  enables  us 
to  obtain  efficient  reduced  models,  i.e.  models  for  which  solution  cost  is  independent 
of  the  full  model  size. 


The  design  problem  of  interest  here  is  to  ford  the  best  materials  and  geometry 
lengths  combination  to  maximize  the  cooling  efficiency.  Note  that  the  output  is 
directly  related  to  the  cooling  efficiency,  that  is,  the  smaller  the  output,  the  better 
the  cooling  efficiency.  Therefore  the  design  problem  becomes  finding  the  best  set  of 
parameters  z  to  minimize  the  output  y  defined  in  (4.12). 
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4.3  Numerical  Results 


In  this  chapter,  the  ranges  (bound  constraints)  of  the  parameters  are  as  follows 
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(4.21) 

1 

to 

<  10, 

3  =  - 

..,4, 

(4.22) 

b^j- 1 

<  10, 

3  =  !,  • 

..,4, 

(4.23) 

h  j 

<  2.5, 

3  =  !» • 

..,4, 

(4.24) 

0.1  < 

bn 
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(4.25) 

It  should  be  pointed  out  that  the  parameter  range  is  wide  in  the  sense  that  the  upper 
bound  for  a  parameter  is  two  orders  of  magnitude  larger  than  the  corresponding  lower 
bound.  Thus,  this  presents  a  challenge  for  model  reduction  methods  in  creating  a 
reduced  model  of  moderate  size  that  accurately  captures  the  full  model  behavior  over 
such  a  wide  range  of  parameter  values  in  multidimensional  parametric  input  space. 
We  also  introduce  the  baseline  parameter  set  as 


baseline 

=  0.4, 

*  =  1,2, 

(4.26) 

baseline 

=  0.6, 

i  =  3, 4, 

(4.27) 
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i  =  5,6, 

(4.28) 

baseline 

=  1-2, 

*  =  7,8, 

(4.29) 
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=  1.0, 
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,16, 

(4.30) 
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(4.31) 
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(4.32) 
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(4.33) 
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^baseline 

°Ai 


=  0.25,  i  —  1, ...  ,4,  (4.35) 

^baseline  =  ^  (4.36) 

so  that  if  any  parameter  is  not  allowed  to  vary,  it  will  take  its  corresponding  baseline 
value. 

For  all  problems  considered  in  the  following,  since  it  is  too  computationally  expen¬ 
sive  to  cover  multidimensional  parameter  spaces  with  full  factorial  search,  we  limit 
ourselves  to  define  the  maximum  output  error  to  be  the  maximum  error  between 
the  full  and  the  reduced  model  outputs  over  a  random  set  of  105  parameters  in  the 
parameter  space  under  consideration.  For  the  model-constrained  adaptive  sampling 
method,  we  choose  e  in  Step  2  of  Algorithm  3.1  to  be  e  =  l.Oe  —  14. 

Initial  guesses  for  the  model-constrained  adaptive  sampling  method  are  obtained 
from  logarithmic  random  sampling.  The  smallest  allowable  distance  between  an  ini¬ 
tial  guess  and  all  existing  sample  locations  is  chosen  to  be  mm.i{zlmax  —  zlmin }.  Unless 
otherwise  specified,  the  reduced  basis  is  computed  by  performing  Gram-Schmidt  or- 
thogonalization  on  the  snapshots.  Finally,  the  computation  time  is  measured  on  a 
dual  core  64-bit  personal  computer  with  3.2GHz  Pentium  processor. 


4.3.1  Error  Indicator  Using  the  Euclidean  ^-norm  versus  the 
Hilbert  Dual  Norm 

As  discussed  in  Section  3.1.2,  one  can  use  the  residual  norm  as  the  error  indicator  in 
the  greedy  optimization  problem.  Intuitively,  the  adaptive  sampling  procedure  will 
drive  the  residual  to  zero  as  more  greedy  cycles  are  taken.  For  steady  problems  that 
are  linear  in  the  state  vector,  using  the  residual  in  the  2-norm ,  for  example,  one  can 
show  that 

l|y  —  yr||2  <  ||C||2||A— 1||2||R||2,  (4.37) 

and  therefore  the  true  output  error  is  also  driven  down  to  zero  as  the  residual  norm 
approaches  zero.  The  question  is  now  which  norm  to  use  for  the  residual.  The  two 
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obvious  choices  are  the  residual  in  the  Euclidean  2-norm 


||R||^  =  RtR,  (4.38) 

and  in  the  dual  norm  [61]  (discrete  H~1  —  norm ) 

||R|||>  =  RtA-1R  (4.39) 

where  A  is  a  SPD  matrix  that  is  given  by,  e.g.  using  the  finite  element  method, 

Aij  =  [  V(fli  ■  Vipj  dfl  +  Y]  /  <Pi<Pj  dTlf,  (4.40) 

k=i  Jnf 

where  ifi  denotes  the  usual  hnite  element  hat  functions.  The  theoretical  problem  with 
the  Euclidean  2-norm  is  that  the  residual  using  this  norm  is  a  sum  of  an  increasing 
number  of  terms  as  the  grid  (or  mesh)  is  refined  [116],  and  hence  theoretically  it  will 
converge  to  infinity  as  the  grid  size  approaches  zero.  The  advantage  of  the  Euclidean 
2-norm  is  that  it  is  simple  and  may  be  the  only  choice  for  problems  that  are  not 
from  a  discretization  of  some  set  of  PDEs,  such  as  a  molecular  dynamic  simulation 
problem.  The  residual  in  the  dual  norm,  on  the  other  hand,  is  more  expensive  to  use 
since  it  involves  the  inversion  of  A,  but  as  the  grid  size  approaches  zero  it  will  not 
approach  infinity  due  to  the  presence  of  A-1  as  the  stabilization. 

In  practice,  the  grid  size  is  hnite  and  we  would  like  to  investigate  the  impact 
of  these  norms  on  the  quality  of  the  reduced  model  using  the  adaptive  sampling 
approach.  In  particular,  we  consider  using  these  norms  for  the  residual  in  the  cost 
functional  in  (3.10)  to  find  a  reduced  basis,  and  hence  a  reduced  model.  We  consider 
a  test  case  of  the  thermal  fin  problem  with  five  parameters,  namely  the  thermal 
conductivities  of  the  four  subfins  and  the  Biot  number,  i.e. 

Z  =  {/ci,  /c2,  «3,  «4,  Biot}.  (4.41) 

Figure  4-3  shows  a  coarse  grid  with  1333  nodes,  a  medium  grid  with  4760  nodes 
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and  a  fine  grid  with  17899  nodes  for  the  thermal  fin.  Recall  that  a  finite  element 
model  with  linear  triangular  elements  is  used  for  the  thermal  fin  problem,  and  hence 
the  number  of  nodes  is  the  same  as  the  number  of  unknowns  (because  there  are  no 
Dirichlet  boundary  conditions  in  this  problem). 

Next,  we  run  the  adaptive  sampling  procedure  with  25  greedy  cycles  using  two 
different  greedy  optimization  problems  that  have  the  same  form  as  (3.10)-(3.13)  but 
with  the  residual  in  the  Euclidean  2-norm  and  the  Hilbert  dual  norm  as  the  cost 
functional.  We  have  used  the  same  sequence  of  initial  guesses  with  25  parameter 
points  for  the  two  optimization  problems.  In  Table  4.1  the  maximum  output  error  is 
used  to  compare  the  quality  of  the  resulting  reduced  models  when  the  grid  is  refined. 
As  can  be  seen,  using  the  adaptive  sampling  procedure  for  the  thermal  fin  problem, 
the  difference  between  using  the  residual  in  the  2-norm  and  in  the  dual  norm  is  small. 
In  fact,  for  a  finite  grid  size,  these  norms  are  equivalent  in  the  sense 

ow||R||2  <  ||R|| D  <  cw||R||2,  (4.42) 

where  amin  and  amax  are  the  maximum  and  the  minimum  eigenvalues  of  the  matrix 
A~h  The  similarity  in  performance  is  further  confirmed  in  Figure  4-4  in  which  we 
plot  the  maximum  output  error  versus  the  number  of  reduced  basis  vectors.  Again, 
the  difference  is  so  small  that,  for  a  finite  grid  size,  either  the  residual  in  the  2-norm 
or  in  the  dual  norm  can  be  used  as  the  error  indicator  in  the  adaptive  sampling 
procedure. 

Table  4.1:  Adaptive  sampling  procedure  for  the  thermal  fin  problem  with  five  param¬ 
eters  using  the  residual  in  the  2-norm  and  in  the  dual  norm  as  the  cost  functional. 
The  maximum  output  error  is  shown  for  all  three  grids. 


Grids 

max  ||  y  yr  II 1 

(S  =  M2) 

max  || y  yr  II 1 

(0  =  ||R||d) 

Coarse 

3.5830e-06 

1.5867e-05 

Medium 

8.9126e-06 

1.5575e-05 

Fine 

7.5282e-06 

4.7335e-06 

(a)  Coarse  grid 
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(b)  Medium  grid 


(c)  Fine  grid 


Figure  4-3:  The  thermal  fin  with  a  coarse  grid  with  1333  nodes,  a  medium  grid  with 
4760  nodes  and  a  fine  grid  with  17899  nodes. 


From  now  on  to  the  end  of  this  chapter,  the  fine  grid  is  used  for  all  problems  and 
only  the  residual  in  the  2-norm  is  used  as  the  error  indicator. 

4.3.2  Performance  of  the  Bound-Constrained  Optimization 
Solver 

In  this  section,  we  investigate  the  efficiency  of  the  optimization  solver  for  the  greedy 
optimization  problem  discussed  in  Section  3.1.3.  In  Figure  4-5 (a),  we  show  the  number 
of  Newton  steps  versus  the  number  of  greedy  cycles  for  the  case  of  34  parameters.  As 
can  be  seen,  the  maximum  number  of  Newton  steps  is  40,  which  occurs  at  the  24th 
greedy  cycle.  For  the  other  adaptive  cycles,  the  number  of  Newton  steps  is  less  than 
30.  Figure  4-5 (b)  shows  the  maximum  number  of  Newton  steps  over  2d  greedy  cycles 
for  cases  with  a  varying  number  of  parameters  d.  It  can  be  seen  that  for  these  cases, 
the  number  of  Newton  steps  is  O(d). 

Recall  that,  due  to  the  combination  of  the  trust  region  strategy  and  the  inexact 
Newton-CG  method,  quadratic  convergence  is  observed  for  all  greedy  cycles  that 
satisfy  the  quadratic  convergence  condition  (i.e.  when  the  iterate  is  close  enough  to  a 
local  maximizer).  Typical  quadratic  convergence  is  shown  in  Table  4.2  in  which  the 
cost  functional  and  the  scaled  gradient  are  shown  along  with  the  number  of  Newton 


(a)  Medium  Grid  (b)  Fine  Grid 


Figure  4-4:  Maximum  reduced  model  output  error  using  the  residual  in  the  2-norm 
versus  in  the  dual  norm  for  the  thermal  problem  with  five  parameters. 


(a)  Number  of  Newton  steps  versus  number  (b)  Maximum  number  of  Newton  steps  ver- 
of  greedy  cycles  for  34-parameter  case  sus  number  of  parameters 


Figure  4-5:  Performance  of  the  bound-constrained  optimization  solver:  a)  number 
of  Newton  steps  versus  number  of  greedy  cycles  for  34-parameter  case;  b)  maximum 
number  of  Newton  steps  versus  number  of  parameters 

steps.  When  the  iterate  is  far  way  from  the  local  maximizer,  the  trust  region  method 
picks  the  steepest  ascent  direction  and  the  step  is  allowed  to  be  as  large  as  possible. 
This  happens  for  the  first  16  Newton  steps.  Once  the  iterate  is  close  to  the  local 
maximizer,  the  Newton  direction  is  picked  and  quadratic  convergence,  e.g.  clearly  for 
the  scaled  gradient  in  Table  4.2,  is  observed  for  the  last  three  Newton  steps.  The  fast 
convergence  can  also  be  seen  in  the  cost  functional  Q  which  increases  by  9  orders  of 
magnitude  from  the  initial  guess  to  the  local  maximizer  within  18  Newton  steps. 

Table  4.2:  Typical  quadratic  convergence  in  the  scaled  gradient  || D 1 1 2  from  the 
bound  constrained  optimization  solver  for  the  case  with  34  parameters.  The  data  are 
shown  for  the  third  greedy  cycle. 


Number  of  Newton  steps 

G 

l|DV£||2 

1 

3.5306e-07 

1.1503e-05 

13 

3.4580e+02 

2.5811e-01 

14 

3.4598e+02 

6.0533e-02 

15 

3.4603e+02 

6.2677e-02 

16 

3.4603e+02 

1.5441e-01 

17 

3.4603e+02 

7.5716e-04 

18 

3.4603e+02 

2.0260e-10 
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4.3.3  Output  Error  Indicator  versus  True  Output  Error 

In  this  section,  we  will  compare  the  offline  cost  as  well  as  the  quality  of  the  reduced 
model  using  the  model-constrained  adaptive  sampling  approach  in  Algorithm  3.1  with 
an  error  indicator  and  the  true  output  error.  For  the  error  indicator,  since  we  only 
compute  the  full  solution  snapshots  at  the  optimal  parameter  points  found  from  the 
greedy  optimization  problem,  only  these  snapshots  are  used  to  enrich  the  reduced 
basis.  For  the  true  error,  we  also  take  only  the  snapshots  at  the  optimal  parameter 
points  to  enrich  the  reduced  basis.  For  both  approaches,  the  same  sequence  of  initial 
guesses  is  used  in  each  greedy  optimization  cycle,  and  Gram-Schmidt  orthogonaliza- 
tion  is  used  to  update  the  reduced  basis. 

The  examples  considered  here  have  11  and  21  parameters,  that  is, 

z  =  {/ci,  ...,ks,  Bi,  b !,...,  65},  (4.43) 

or 

z  =  {/Cl, . . . ,  Kio,  Bi,  610}.  (4.44) 

We  use  the  number  of  full  matrix  factorizations  as  the  measure  for  the  offline  cost 
to  compare  the  quality  of  the  reduced  models,  since  this  is  the  dominant  cost  of  the 
reduction  algorithm.  For  the  11-parameter  case,  Figure  4-6(a)  shows  that  the  required 
number  of  matrix  factorizations  to  reach  a  given  error  level  is  approximately  an  order 
of  magnitude  larger  for  the  true-error  approach;  however,  for  the  same  number  of 
basis  functions  retained  in  the  reduced  basis,  Figure  4-6 (b)  shows  that  using  the  true 
error  rather  than  the  indicator  leads  to  more  efficient  (i.e.  smaller  for  a  given  error 
level)  reduced  models.  This  result  might  be  intuitively  expected,  since  the  optimal 
parameter  points  based  on  the  true  output  error  should  better  target  reduction  of 
the  true  output  error  than  those  points  selected  using  the  error  indicator.  However, 
there  is  no  guarantee  that  this  will  always  be  the  case,  as  shown  in  Figure  4-7  for  the 
case  of  21  parameters.  Figure  4-7(a)  shows  that  the  number  of  matrix  factorizations 
is  again  about  an  order  of  magnitude  larger  for  the  true-error  approach.  For  smaller 
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basis  sizes,  Figure  4-7(b)  shows  that  the  output  error  is  again  smaller  than  for  models 
obtained  using  the  error  indicator;  however,  for  larger  basis  sizes,  the  true-error  and 
error-indicator  approaches  give  equivalently  good  reduced  models. 

Figures  4-6  and  4-7  demonstrate  a  general  tradeoff  in  the  model-constrained  sam¬ 
pling  methodology:  if  one  is  willing  to  invest  larger  offline  cost  to  compute  the  reduced 
model,  then  using  the  true  error  to  select  the  parameter  sample  points  can  lead  to 
more  efficient  models.  For  some  problems,  such  as  real-time  applications,  minimiz¬ 
ing  the  size  of  the  reduced  model  may  be  critical;  in  that  case,  one  might  choose  to 
use  the  true-error  approach.  For  very  large-scale  problems,  however,  the  cost  of  the 
true-error  approach  may  be  prohibitively  high;  in  that  case,  the  error  indicator  is  an 
attractive  option.  In  many  of  the  numerical  experiments  performed  for  the  thermal 
fin  problem,  and  as  demonstrated  by  the  results  in  Figure  4-7,  the  difference  in  quality 
between  the  true-error  and  error-indicator  sample  points  tends  to  be  larger  in  early 
greedy  cycles.  Since  the  error  function  becomes  more  multimodal  as  the  number  of 
greedy  cycles  increases,  the  chances  of  sampling  at  a  local  maximum  are  increased, 
and  thus  the  difference  between  the  error- indicator  and  true-error  approaches  may 
not  be  as  great.  One  could  therefore  also  conceive  of  using  a  combination  of  the  two 
error  metrics,  i.e.  using  the  true  error  for  early  greedy  cycles  and  the  error  indicator 
for  later  cycles,  in  an  attempt  to  balance  offline  cost  with  the  quality  of  the  reduced 
model. 

As  discussed  in  Section  3.1.6,  when  using  the  true  output  error,  one  has  interme¬ 
diate  full  state  snapshots  and  full  adjoint  snapshots  available  at  all  Newton  steps,  i.e. 
the  available  snapshot  information  includes  but  is  not  limited  to  those  solutions  at  the 
optimal  parameter  points.  We  compare  four  methods  of  updating  the  reduced  basis, 
described  in  Section  3.1.6  and  summarized  in  Table  4.3.  In  the  first  method,  we  add 
all  the  snapshots  of  the  current  greedy  cycle  to  the  current  reduced  basis  using  the 
Gram-Schmidt  procedure  (i.e.  adding  only  information  that  is  linearly  independent  of 
the  current  reduced  basis  vectors).  In  the  second  method,  we  store  all  the  snapshots, 
which  is  possible  for  the  thermal  fin  problem  for  a  small  number  of  greedy  cycles,  and 
then  perform  the  POD  method  on  the  complete  snapshot  set  (i.e.  that  comprising 
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(a)  Maximum  error  versus  the  number  of  full  matrix  factorizations 


(b)  Maximum  error  versus  the  number  of  reduced  basis  vectors 


Figure  4-6:  Error  indicator  versus  true  output  error  for  the  thermal  fin  with  11 
parameters.  The  same  sequence  of  initial  guesses  is  used  for  both  true-error  and 
error-indicator  approaches,  and  the  Gram-Schmidt  procedure  is  used  to  update  the 
reduced  basis. 
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(a)  Maximum  error  versus  the  number  of  full  matrix  factorizations 


(b)  Maximum  error  versus  the  number  of  reduced  basis  vectors 


Figure  4-7:  Error  indicator  versus  true  output  error  for  the  thermal  fin  with  21 
parameters.  The  same  sequence  of  initial  guesses  are  used  for  both  true-error  and 
error-indicator  approaches,  and  the  Gram-Schmidt  procedure  is  used  to  update  the 
reduced  basis. 
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the  newly  computed  snapshots  of  the  current  greedy  cycle  and  all  the  snapshots  from 
the  previous  greedy  cycles).  Denote  x%  to  be  the  singular  value  corresponding  to  the 
jth  pop  pasis  vector,  which  is  added  to  the  reduced  basis  if  the  following  condition 
is  satisfied  [40,51] 

—  <  V,  (4.45) 

1  Xj 

where  ns  is  the  number  of  snapshots,  and  rj  <  1.  We  choose  rj  =  0.99999999  for 
the  second  method.  The  third  method  is  to  perform  the  POD  method  on  the  newly 
computed  snapshots  of  the  current  greedy  cycle  with  rj  =  0.99999999,  and  then  add 
only  the  dominant  POD  basis  vectors  to  the  current  reduced  basis  using  the  Gram- 
Schmidt  procedure.  The  fourth  method  considered  is  the  base  case  employed  in 
Figures  4-6  and  4-7,  i.e.  using  the  Gram-Schmidt  procedure  only  on  the  snapshots  at 
the  optimal  parameter  points. 


Table  4.3:  Four  methods  of  updating  the  reduced  basis  when  using  the  true  output 
error. 


Method 

Intermediate 
information  included 

Basis  updating 

TrueErrorGS 

Yes 

Gram-Schmidt 

TrueErrorPOD 

Yes 

POD 

TrueErrorPOD-GS 

Yes 

POD  and  then  Gram-Schmidt 

TrueErrorGS-Opt 

No 

Gram-Schmidt 

For  these  four  methods,  the  reduced  bases  at  each  greedy  cycle  are  not  the  same, 
and  hence  the  snapshots  found  in  each  cycle  are  different.  Figures  4-8 (a)  and  4-9 (a) 
show  that,  with  a  same  number  of  matrix  factorizations,  the  first  three  methods  yield 
reduced  models  of  equivalent  accuracy.  However,  as  can  be  seen  in  Figures  4-8  (b) 
and  4-9  (b),  with  the  same  number  of  reduced  basis  vectors,  the  first  method  leads 
to  reduced  models  that  are  less  accurate.  That  is,  adding  all  new  snapshots  leads 
to  reduced  models  that  are  inefficient  and  unnecessarily  large.  The  comparison  with 
the  base  approach  shows  that  discarding  the  intermediate  information  leads  to  higher 
offline  cost,  but  significantly  smaller  reduced  models  for  a  given  level  of  output  error. 
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(a)  Maximum  error  versus  the  number  of  full  matrix  factorizations 


(b)  Maximum  error  versus  the  number  of  reduced  basis  vectors 


Figure  4-8:  Different  methods  of  updating  the  reduced  basis  for  the  case  with  11 
parameters.  The  comparisons  are  the  true  output  error  using  the  Gram-Schmidt 
procedure  for  all  snapshots  (  TrueErrorGS),  the  true  output  error  with  the  POD 
method  for  all  snapshots  after  each  greedy  cycle  (TrueErrorPOD),  the  true  output 
error  with  the  Gram-Schmidt  procedure  on  the  POD  vectors  of  the  snapshots  of  the 
current  greedy  cycle  and  on  the  current  reduced  basis  (TrueErrorPOD-GS),  and  the 
true  output  error  with  the  Gram-Schmidt  procedure  only  on  the  optimal  parameter 
points  (TrueErrorGS-Opt). 
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(a)  Maximum  error  versus  the  number  of  full  matrix  factorizations 


(b)  Maximum  error  versus  the  number  of  reduced  basis  vectors 


Figure  4-9:  Different  methods  of  updating  the  reduced  basis  for  the  case  with  21 
parameters.  The  comparisons  are  the  true  output  error  using  the  Gram-Schmidt 
procedure  for  all  snapshots  (  TrueErrorGS),  the  true  output  error  with  the  POD 
method  for  all  snapshots  after  each  greedy  cycle  (TrueErrorPOD),  the  true  output 
error  with  the  Gram-Schmidt  procedure  on  the  POD  vectors  of  the  snapshots  of  the 
current  greedy  cycle  and  on  the  current  reduced  basis  (TrueErrorPOD-GS),  and  the 
true  output  error  with  the  Gram-Schmidt  procedure  only  on  the  optimal  parameter 
points  (TrueErrorGS-Opt). 
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In  summary,  either  the  error  indicator  or  the  true  error  can  be  used  as  the  cost 
functional  in  the  model-constrained  adaptive  sampling  approach.  For  the  thermal 
fin  problem,  the  error- indicator  approach  leads  to  accurate  reduced  models  with  ap¬ 
proximately  one  order  of  magnitude  reduction  in  offline  cost  compared  to  using  the 
true  error.  Unless  it  is  critical  to  decrease  the  size  of  the  reduced  model  as  much  as 
possible,  using  the  error  indicator  is  the  recommended  approach,  especially  for  very 
large-scale  systems  where  the  additional  offline  computational  cost  of  using  the  true 
error  may  be  prohibitively  high.  Further,  if  using  the  true-error  approach,  although 
intermediate  state  and  adjoint  information  is  available  and  can  be  included  in  the 
basis  updating  process,  the  results  indicate  that  doing  so  compromises  the  quality  of 
the  resulting  reduced  models. 

In  the  next  section,  we  compare  the  model-constrained  adaptive  sampling  method 
with  other  sampling  approaches.  For  all  results  that  follow,  the  error- indicator  ap¬ 
proach  is  used. 

4.3.4  Model-Constrained  Adaptive  Sampling  versus  Other 
Sampling  Methods 

In  this  section  we  compare  the  model-constrained  adaptive  sampling  approach  using 
the  residual  error  indicator  as  the  objective  function  with  other  sampling  approaches. 
We  first  consider  the  greedy  sampling  method,  which  is  the  fundamental  concept 
underlying  the  model-constrained  approach.  We  then  compare  the  model-constrained 
approach  with  four  statistically-based  sampling  methods. 

Model-Constrained  Adaptive  Sampling  Approach  versus  the  Greedy  Sam¬ 
pling  Approach 

To  apply  the  greedy  sampling  method  in  [61,62,68],  one  needs  to  determine  a  train¬ 
ing  parameter  set  with  ntram  parameters.  At  each  point  in  this  parameter  set,  the 
reduced  states  are  computed  by  solving  the  reduced  model.  The  error  estimator,  i.e. 
here  the  residual  in  the  2-norm,  is  then  computed.  The  location  in  the  training  pa- 
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rameter  set  at  which  the  error  estimator  is  maximal  is  then  found,  the  reduced  basis 
is  updated  with  the  snapshot  at  this  sample  location,  and  the  updated  reduced  model 
is  computed.  To  generate  the  training  set  for  the  greedy  sampling  approach,  we  use 
two  different  methods,  namely  logarithmic  random  sampling  and  Latin  hypercube 
sampling  (LHS). 

For  the  model-constrained  adaptive  sampling  approach,  we  use  logarithmic  ran¬ 
dom  sampling  to  generate  initial  guesses  for  the  greedy  optimization  problem.  The 
residual  in  the  2-norm  is  used  as  the  error  indicator;  hence,  only  snapshots  at  the  op¬ 
timal  parameter  points  are  computed.  A  total  of  100  greedy  cycles  are  used  for  both 
model-constrained  sampling  and  greedy  sampling  methods;  hence,  both  methods  will 
provide  100  sample  points  at  which  the  snapshots  are  computed  to  form  the  reduced 
basis.  This  comparison  set-up  gives  us  the  same  offline  cost  corresponding  to  the 
full  model  solve,  costF,  in  (3.25)  for  both  methods.  The  only  difference  is  the  offline 
cost  corresponding  to  the  reduced  model  solve,  cost^.  For  the  model-constrained 
approach,  cost^  for  steady  problems  is  approximately  given  as 

n°  2 

costlde,constrained  ~  d  Y  o  f  +  [2(x  +  nHv)  +  n%]j 2,  (4.46) 

3= 1  ^ 

where,  following  the  affine  decomposition  in  (4.14),  we  have  also  incorporated  the 
cost  of  forming  the  reduced  matrix  Ar  as  n|,i2,  where  n©  =  48  is  the  number  of  terms 
in  (4.14).  The  offline  cost  corresponding  to  the  reduced  model  solve  for  the  greedy 
sampling  approach  (either  with  logarithmic  random  sampling  or  LHS)  is  given  by 

n°  2 

cost  Jeedy  ~  ntrain  Y  3  f  +  (ne  +  l)f-  (4-47) 

3= 1 

At  each  greedy  cycle,  the  greedy  sampling  approach  needs  to  solve  ntrain  dense  reduced 
models  at  ntrain  parameters  to  find  the  reduced  states  in  order  to  compute  the  error 
estimator.  The  model-constrained  approach  needs  to  solve  the  greedy  optimization 
problem  in  each  greedy  cycle.  The  difference  in  the  cost  in  each  greedy  cycle  is 
therefore  the  multiplicative  constant  d  for  the  model-constrained  approach  and  ntra in 


for  greedy  sampling  methods.  For  the  following  comparisons,  we  choose  ntram  = 
104  >  d  E  {11,21}. 

Figures  4-10(a)  and  4-10(b)  show  the  comparison  for  the  three  methods  for  the 
case  of  11  and  21  parameters,  respectively.  It  can  be  seen  that  the  maximum  output 
error  obtained  using  the  model-constrained  approach  is,  for  the  most  part,  comparable 
to  that  of  the  greedy  sampling  method  with  logarithmic  random  training  points.  For 
the  case  with  21  parameters,  the  model-constrained  method  is  able  to  achieve  an 
order  of  magnitude  better  error  than  the  greedy  approach  for  larger  reduced  model 
sizes.  Using  the  greedy  sampling  method  with  LHS  training  points  led  to  larger  errors 
than  those  obtained  using  the  other  two  methods. 

A  key  difference  between  the  methods  is  that  the  greedy  sampling  approach  finds 
globally-optimal  parameter  points  within  the  discrete  training  set  (via  exhaustive 
search)  while  the  model-constrained  approach  finds  locally-optimal  parameter  points 
in  the  continuous  parameter  space  (by  solving  an  optimization  problem).  As  a  result, 
unless  ntram  is  sufficiently  large,  the  training  set  may  not  adequately  cover  important 
parameter  regions,  particularly  as  the  dimension  of  the  parametric  space  becomes 
large.  This  difference  is  highlighted  by  the  results  in  Figure  4-10(b),  where  even  104 
training  points  were  not  sufficient  for  the  greedy  sampling  method  to  find  near-optimal 
sample  points  in  later  greedy  iterations. 

Although  the  total  number  of  large-scale  matrix  factorizations  is  100  for  both 
the  model-constrained  and  greedy  sampling  methods,  the  actual  offline  cost  differs 
substantially  between  the  two.  Table  4.4  compares  the  CPU  time  required  to  compute 
the  reduced  basis  for  the  three  approaches  for  the  case  of  21  parameters.  It  can  be 
seen  that  the  model-constrained  approach  is  approximately  16  times  faster  than  the 
greedy  sampling  approaches.  This  difference  is  due  to  the  need  for  exhaustive  search 
over  104  training  points  on  every  greedy  iteration  in  the  greedy  sampling  method, 
which  for  this  case  is  the  dominant  offline  cost.  This  could  be  a  particular  concern 
as  the  number  of  parameters,  d,  and  hence  the  necessary  number  of  training  points, 
Strain,  increases. 
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(b)  21-parameter  case 


Figure  4-10:  Model-constrained  sampling  approach  versus  greedy  sampling  ap¬ 
proaches  (Greedy-LogRandom:  greedy  sampling  with  logarithmic  random  training 
parameter  set,  and  Greedy-LHS:  greedy  sampling  with  LHS  training  parameter  set) 
over  100  greedy  cycles.  A  training  set  with  104  training  points  is  used  in  the  greedy 
search  for  the  greedy  sampling  approaches. 
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Table  4.4:  The  offline  cost  in  CPU  time  of  the  mo  del- constrained  sampling  approach 
and  the  greedy  sampling  approaches  for  the  case  of  21  parameters.  100  greedy  cycles 
are  taken  for  all  methods  and  ntram  =  104  for  the  greedy  sampling  approaches. 


Model-constrained 

Greedy  sampling 

Greedy  sampling 

sampling 

with  logarithmic  random 

with  LHS 

CPU  time 

0.58  hours 

8.33  hours 

8.33  hours 

Model-Constrained  Adaptive  Sampling  versus  Other  Sampling  Methods 

Next,  we  compare  the  model-constrained  sampling  method  with  statistically-based 
sampling  methods  in  the  context  of  snapshot  generation  for  model  reduction.  In 
particular,  we  compare  our  model-constrained  sampling  with  LHS  sampling,  uniform 
random  sampling,  logarithmic  random  sampling,  and  CVT  sampling.  For  all  methods, 
we  take  100  sample  points  to  generate  100  snapshots,  which  are  then  used  to  form  the 
reduced  basis.  As  can  be  seen  in  Figures  4-11  (a)  and  (b)  for  11  and  21  parameters, 
respectively,  the  model-constrained  sampling  method  outperforms  the  other  methods 
in  the  sense  that,  for  a  given  basis  size,  the  reduced  model  error  is  several  orders  of 
magnitude  smaller  than  that  obtained  using  the  other  methods.  Furthermore,  going 
from  11  to  21  parameters,  the  difference  in  accuracy  of  the  reduced  model  using  model- 
constrained  sample  points  and  those  of  other  methods  is  larger.  This  reflects  the  fact 
that  the  model-constrained  sampling  is  a  model-based  sampling  method;  that  is,  the 
parameter  is  sampled  where  the  indicator  of  the  error  between  the  full  and  the  reduced 
models  is  locally  largest.  The  other  methods,  on  the  other  hand,  use  no  knowledge 
of  the  underlying  model  when  selecting  their  sample  points.  As  the  dimension  of  the 
parametric  space  increases,  it  becomes  more  and  more  difficult  to  adequately  cover 
the  space  with  a  reasonable  number  of  samples  using  the  statistically-based  sampling 
methods. 

Generating  the  sample  points  using  either  logarithmic  random  sampling  or  uniform 
random  sampling  is  very  rapid.  For  the  other  methods,  there  is  additional  overhead 
in  LHS  (due  to  stratification  and  grouping  processes),  in  CVT  sampling  (due  to  its 
iterative  nature),  and  in  model-constrained  sampling  (due  to  solving  the  optimization 
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problem) .  We  also  note  that  while  logarithmic  random  sampling  is  less  expensive  (in 
terms  of  overhead)  than  the  CVT  and  LHS  sampling  methods,  it  leads  to  more 
accurate  reduced  models  in  the  case  of  the  thermal  hn  problem. 

For  the  statistically-based  methods  that  have  some  random  elements  (i.e.  LHS, 
uniform  random  sampling  and  logarithmic  random  sampling),  the  comparison  in  Fig¬ 
ure  4-11  is  limited  to  one  random  set  of  parameters.  A  different  random  draw  could 
lead  to  different  results  and  different  reduced  model  performance.  However,  first  we 
note  that  the  relative  performance  difference  of  the  mo  del- constrained  sampling  ap¬ 
proach  is  significant;  thus,  variation  in  performance  due  to  variability  in  the  random 
draw  is  not  expected  to  alter  the  conclusion  that  the  model-constrained  method  per¬ 
forms  better  than  the  other  methods.  Second,  even  if  one  particular  random  draw 
could  lead  to  a  reduced  model  that  outperforms  the  model-constrained  approach,  it 
is  preferable  to  have  a  systematic  methodology  with  repeatable  performance,  rather 
than  relying  on  a  favorable  random  sample  to  achieve  good  results. 

Next,  we  study  the  sensitivity  of  the  quality  of  model-constrained  sample  points 
using  different  methods  to  generate  the  initial  guesses.  In  particular,  we  take  the 
points  in  parameter  space  corresponding  to  the  logarithmic  random,  uniform  ran¬ 
dom,  LHS,  and  CVT  samplings  that  were  previously  used  to  obtain  the  results  in 
Figure  4-11  as  initial  guesses,  z°,  for  the  model-constrained  sampling  method.  For 
each  parameter  point,  we  then  solve  an  optimization  problem  to  determine  the  cor¬ 
responding  z*,  which  becomes  the  new  sample  point.  Since  our  algorithm  uses  the 
error  indicator,  we  only  perform  a  large-scale  solve  at  the  optimal  parameter  point  z*; 
hence  the  cost  is  not  increased  (with  the  exception  of  the  optimization  solver  over¬ 
head).  Figures  4- 12  (a)  and  (b)  show  the  maximum  output  error  versus  the  number 
of  reduced  basis  vectors  for  both  the  11-  and  21-parameter  cases.  It  can  be  seen  that 
the  quality  of  the  reduced  model  resulting  from  the  model-constrained  approach  is 
relatively  insensitive  to  the  initial  guesses,  at  least  for  the  thermal  hn  problem  con¬ 
sidered  here.  That  is,  the  quality  of  the  sample  points  using  different  methods  to 
generate  initial  guesses  is  more  or  less  the  same.  Again,  this  emphasizes  the  benefit 
of  a  systematic,  repeatable  point  selection  criterion. 
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(b)  21-parameter  case 


Figure  4-11:  Model-constrained  sampling  versus  LHS,  uniform  random  sampling 
(URandom),  logarithmic  random  sampling  (LogRandom)  and  CVT  sampling.  100 
sample  points  are  used  for  all  methods. 
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(a)  11-parameter  case 


(b)  21-parameter  case 


Figure  4-12:  Model-constrained  sampling  method  using  logarithmic  random  sam¬ 
pling  (ModelConstrained-LogRandom),  LHS  (ModelConstrained-LHS),  uniform  ran¬ 
dom  (ModelConstrained-URandom),  and  CVT  (ModelConstrained-CVT)  to  generate 
initial  parameter  guesses,  z°.  100  sample  points  are  used  for  all  methods  of  initial¬ 
ization. 
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4.3.5  Optimal  Design  Application 


In  this  section,  the  reduced  model  is  used  as  a  surrogate  model  in  the  optimal  design 
context.  That  is,  instead  of  performing  the  design  task  on  the  full  model,  the  design 
task  is  carried  out  on  the  reduced  model,  which  is  much  less  expensive.  Recall  that 
the  design  problem  is  to  find  a  combination  of  thermal  conductivities,  Biot  number, 
and  lengths  to  minimize  the  output  average  temperature  (and  hence  maximize  the 
cooling  efficiency). 

Here  we  consider  the  thermal  ffii  optimal  design  problem  with  11  parameters,  as 
shown  in  Table  4.5.  The  reduced  model  is  generated  using  the  model-constrained 
approach  as  in  Section  4.3.4  with  100  basis  vectors.  Using  the  same  initial  design, 
which  is  shown  in  the  second  column  in  Table  4.5,  we  obtain  different  optimal  solutions 
using  the  full  and  the  reduced  model,  as  shown  in  the  third  and  fourth  columns  of 
Table  4.5,  respectively.  Although  the  values  of  the  objective  function  for  these  two 
designs  are  similar,  the  optimal  solution  obtained  from  the  reduced  model  does  not 
satisfy  the  optimality  conditions  of  the  full  problem. 

To  improve  the  quality  of  the  reduced  model  we  run  50  additional  greedy  cycles 
to  obtain  50  more  basis  vectors.  The  optimal  design  for  the  reduced  model  with  150 
basis  vectors  is  shown  in  the  fifth  column  of  Table  4.5.  Again,  the  optimal  solution 
of  the  reduced  model  is  different  from  that  of  the  full  one,  although  the  objective 
function  value  is  the  same.  In  this  case,  the  optimum  found  by  the  reduced  model  is 
a  local  minimum  of  the  full  model. 

The  reduction  factor  in  model  size  is  119  (from  the  full  model  with  17899  states 
to  the  reduced  model  with  150  states),  but  the  speed-up  factor  in  solving  the  optimal 
design  problem  is  58.12/2.49  ~  23  as  shown  in  Table  4.5.  This  is  because  the  full 
system  matrix  is  sparse  (A  €  jpp 7899x17899  }ias  120,603  nonzero  entries)  while  the  re¬ 
duced  matrix  is  not  (Ar  G  1R150x150  has  22,500  nonzero  entries).  For  this  problem  with 
17899  states  (which  is  relatively  small),  the  offline  cost  of  determining  the  reduced 
model  is  not  offset  by  computational  savings  in  solving  the  optimal  design  problem, 
unless  there  is  a  need  to  solve  the  optimization  problem  in  real  time  (in  which  case 
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the  reduced  model  is  essential).  Chapter  6  presents  a  large-scale  probabilistic  analysis 
application.  It  will  be  shown  that  in  that  case,  the  trade-off  between  the  online  and 
offline  cost  is  much  more  dramatic. 

Table  4.5:  Optimal  design  for  the  thermal  fin  problem  with  11  parameters:  full  model 
(17899  states)  versus  reduced  models  with  100  and  150  states. 


Parameters 

Initial  design 

Optimal  solution 
using  the 
full  model 

Optimal  solution 
with  100 
basis  vectors 

Optimal  solution 
with  150 
basis  vectors 

K 1 

2.8657 

10.0 

0.10 

10.0 

2.1151 

10.0 

10.0 

10.0 

0.1888 

10.0 

0.1 

0.1 

K4 

0.2349 

5.25 

10.0 

10.0 

k5 

0.1947 

5.25 

10.0 

0.1 

Bi 

2.5531 

5.0 

5.0 

5.0 

h 

0.2529 

5.0 

5.0 

5.0 

b2 

0.3324 

10.0 

10.0 

10.0 

h 

8.2753 

10.0 

10.0 

10.0 

h 

0.7236 

2.5 

2.5 

2.5 

b>5 

1.6142 

5.0 

5.0 

5.0 

Optimal 

output 

2.678e-03 

2.673e-03 

2.678e-03 

CPU  time 

58.12s 

6.15s 

2.49s 
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Chapter  5 


Computational  Fluid  Dynamic 
Models 


Recent  developments  in  the  field  of  CFD  have  led  to  the  use  of  higher-order  finite 
element  discretizations  for  PDEs.  These  schemes  have  advantages  over  traditional 
finite-volume  methods  by  introducing  higher-order  accuracy  compactly  within  grid 
elements  and  thus  providing  a  significant  decrease  in  the  computational  cost  to  obtain 
reliably  accurate  solutions.  A  Discontinuous  Galerkin  (DG)  formulation  is  used  in  this 
work.  The  unsteady  flow  solver  described  in  this  thesis  is  part  of  a  larger  effort  that 
includes  an  adaptive  meshing  utility,  a  multigrid  solution  algorithm,  gradient-based 
optimization  capability,  and  high-order  visualization  [117]. 

In  this  chapter,  a  steady  DG  formulation  for  the  Euler  equations  is  first  outlined. 
A  linearized  unsteady  DG  formulation  is  then  presented,  and  validated  via  numerical 
comparisons  with  experimental  data.  Next,  a  linearized  model  for  incorporating  geo¬ 
metric  variability  into  the  unsteady  CFD  model  is  described.  Finally,  the  applicability 
of  the  linearized  geometric  variability  model  is  assessed  via  numerical  experiments. 

5.1  Steady  CFD  Model 

In  this  section  we  briefly  review  the  DG  discretization  and  the  solution  method  for 
the  two-dimensional  Euler  equations  as  in  Cockburn  and  Shu  [118]  and  Fidkowski 
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and  Darmofal  [119].  The  two-dimensional  Euler  equations  are  given  by: 


<9w 

~dt 


V  •  JF(w)  =  0, 


where  w  is  the  conservative  state  vector, 


w  = 


t  p  ' 

pu 
pv 

\pEJ 


and  T  =  (Fx,  Fy)  is  the  inviscid  Euler  flux 


= 


(  pu  ^ 
pu 2  +  P 
puv 
\  P^H  ) 


( 


p  .y  _ 


pv 
puv 

pv2  +  P 

V  PVH  J 


(5.1) 


(5.2) 


(5.3) 


In  the  above  equations,  p  is  the  density,  u  and  v  are  respectively  the  x—  and  y— component 
of  velocity,  E  is  the  energy,  P  is  the  pressure,  and  H  =  E  +  P/p  is  the  total  enthalpy. 
The  equation  of  state  is 


P  =  (7  - 1) 


pE  -  h  („*  +  „*) 


(5.4) 


where  7  is  the  ratio  of  specific  heats. 


As  in  the  continuous  finite  element  method,  the  first  step  in  the  DG  method  is 
to  discretize  the  domain  under  consideration,  12,  into  elements  12e.  Next,  a  space  of 
polynomials  of  degree  at  most  p,  W^(12e),  is  defined  on  each  element,  where  h  denotes  a 
representative  element  size  for  the  discretization  (e.g.  the  size  of  the  smallest  element). 
On  each  element  12e,  the  approximate  solution  w/-(  can  be  found  by  enforcing  the 


108 


nonlinear  conservation  law  (5.1)  locally,  for  all  test  functions 

+  [  K)T^(w+  Wfc,n)ds 

+  [  (vft)r^bd(w+,w^,n)cis  =  0,  (5.5) 

Jaoenan 

where  <912  and  <912e  are  the  boundaries  of  the  entire  domain  hi  and  the  element  12e, 
respectively,  and  n  denotes  the  outward-pointing  normal  on  the  boundaries  of  the 
element.  The  terms  w^,  n)  and  7^bd  (W/T, w^,  n)  are  numerical  flux  functions 

for  interior  and  boundary  edges,  respectively,  where  ()+  and  ()“  denote  values  taken 
from  the  interior  and  exterior  of  the  element.  The  interior  flux  function  is  computed 
using  the  Roe-averaged  flux  function  [120]  and  contributes  over  element  boundaries 
that  do  not  belong  to  the  domain  boundary,  denoted  by  <912e\<91 1  The  fluxes  on  the 
common  boundaries  of  <91 2e  and  <911,  denoted  by  <9 12e  fl  <912,  are  computed  using  the 
inner  state  and  boundary  condition  data. 

The  final  form  of  the  DG  discretization  is  constructed  by  selecting  a  basis  for 
W^(12e).  The  approximate  solution  wh  on  each  element  is  assumed  to  be  a  linear 
combination  of  the  basis  functions  4>j, 

nb 

w  h(t,x,y)  =  ■j(t)(/>j(x,y)^  (5.6) 

j= i 

where  w j(t)  gives  the  modal  content  of  4>j  on  element  12e,  and  rib  is  the  number  of 
basis  functions  required  to  describe  12e)  (e.g.  rib  =  1  for  p  =  0  and  rib  =  3  for 
p  =  1).  The  complete  set  of  unknown  quantities  for  the  DG  formulation  comprises 
the  values  of  w j(t)  for  every  element  in  the  spatial  domain.  These  quantities  will  be 
contained  in  the  vector  w  £  1R™,  where  n  is  the  total  number  of  unknowns,  which 
depends  both  on  the  number  of  elements  in  the  discretization  and  on  the  polynomial 
order  p. 

For  steady-state  flows,  pseudo  time-stepping  is  used  to  improve  the  initial  transient 
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behavior  of  the  solver.  A  backward  Euler  discretization  in  time  is  used  so  that  the 
final  discrete  equations  are 

EXt  (™n+1  _  ™n)  +  Q(wn+1)  = 0  (5.7) 

where  At  is  the  timestep,  wn  is  the  solution  w  at  a  time  tn,  E  is  the  mass  matrix, 
and  Q  is  the  vector  representing  the  final  three  terms  of  (5.5).  This  nonlinear  system 
is  solved  using  a  p-multigrid  scheme  with  a  line  Jacobi  smoother  [117,119]. 


5.2  Unsteady  CFD  Model 

The  unsteady  Euler  equations  using  the  DG  spatial  discretization  can  be  written 

E^  +  Q(w,u)  =  °  (5.8) 


where  u(t)  G  1RP  is  a  vector  containing  p  external  forcing  inputs,  such  as  prescribed 
motion  of  the  domain  boundary  or  incoming  flow  disturbances.  In  addition,  we  define 
a  set  of  q  output  quantities  of  interest,  contained  in  the  vector  y  G  !Rg  and  defined 
by  the  nonlinear  function  V 

y  =  V(w).  (5.9) 


For  unsteady  computations,  a  second-order  backward  Euler  temporal  discretiza¬ 
tion  is  applied  to  (5.8).  The  resulting  nonlinear  equations  are  solved  using  a  Newton 
solver.  Grid  motion  is  implemented  using  a  simple  Jacobi  smoothing  formulation.  The 
motion  of  grid  point  j  is  defined  by  the  change  in  x—  and  y— coordinates,  ( Sxj,5yj ), 
and  computed  as 


6xj 


1 

r 


5xk , 

keNj 


(5.10) 


where  r  is  the  number  of  the  neighbors  chosen  to  influence  the  grid  point  j  and  Nj  is 
the  set  containing  the  corresponding  set  of  r  neighboring  points.  Larger  values  of  r 
lead  to  increased  grid  motion  and  smoother  grids.  The  motion  of  grid  points  on  the 
domain  boundary  is  prescribed  according  to  the  corresponding  external  input  (e.g. 
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prescribed  motion  of  an  airfoil). 


In  many  cases  of  interest,  the  unsteady  flow  solution  can  be  assumed  to  be  a  small 
perturbation  from  steady-state  conditions.  This  allows  the  unsteady  governing  equa¬ 
tions  to  be  linearized,  which  reduces  the  computational  cost  of  solution  considerably. 
The  linearized  version  of  equations  (5.8,  5.9)  can  be  written  in  standard  state-space 
form,  which  is  a  special  case  of  the  general  model  (2.4)  that  is  obtained  from  the  DG 
CFD  linearization, 

(J-V 

E—  =  Ax  +  Bu,  y  =  Cx,  (5.11) 


where  x  G  IR"  is  the  state  vector  containing  the  n  perturbations  in  flow  unknowns 
from  the  steady-state  solution  wss,  that  is  w (t)  =  wss+x(f).  The  matrices  A  G  Rnxn, 
B  G  IRnxp,  and  C  G  !R9Xn  in  (5.11)  have  constant  coefficients  evaluated  at  steady- 
state  conditions  and  arise  from  the  linearization  of  (5.8)  and  (5.9)  as  follows. 


9R  dR  r_^_ 

--fa’ 


(5.12) 


By  considering  harmonic  inputs  at  a  frequency  u,  u  =  u.e^ut,  the  linearized  equa¬ 
tions  (5.11)  can  also  be  written  in  the  frequency  domain  as 


[j'cjE  —  A]  x  =  Bu,  y  =  Cx, 

where  x  =  xe^4  and  y  =  yelut. 


(5.13) 


5.3  CFD  Model  Validation 

Results  are  presented  for  two  unsteady  examples:  a  NACA  0012  airfoil  and  the  first 
standard  cascade  configuration.  Results  are  shown  to  validate  the  unsteady  CFD 
models  by  comparison  with  airfoil  and  cascade  experimental  data. 


Ill 


5.3.1  NACA  0012  Airfoil  Example 


Figure  5-1  shows  nonlinear  and  linearized  CFD  results  compared  with  experimental 
data  [1]  for  a  NACA  0012  airfoil  in  rigid  pitching  motion.  The  steady-state  flow 
has  a  Mach  number  of  0.6  and  angle  of  attack  of  2.89°.  The  unsteady  simulation 
is  then  carried  out  with  the  pitching  input  a(t)  =  2.41  sin(0.8874t),  where  a(t)  is 
the  perturbation  angle  about  the  steady  state  angle  of  attack.  Figure  5-1  shows 
the  pressure  coefficient  distribution  at  one  particular  instant  in  time;  it  can  be  seen 
that  the  nonlinear  CFD  predictions  match  well  with  the  experimental  data.  There 
is  a  region  with  a  weak  shock  on  the  upper  surface  that  the  linearized  code  cannot 
resolve.  Elsewhere  the  agreement  between  the  experimental  data  and  linearized  model 
is  acceptable. 


Figure  5-1:  Comparison  of  CFD  predictions  and  experimental  data  [1]  for  NACA 
0012  airfoil  in  unsteady  pitching  motion.  The  pressure  coefficient  distribution  on  the 
airfoil  surface  is  shown  for  t  =  0.00311. 


5.3.2  The  First  Standard  Cascade  Configuration  Example 

Unsteady  computations  for  cascade  flows  can  be  carried  out  efficiently  by  exploiting 
spatial  periodicity  and  linearity.  By  working  with  the  frequency  domain  equations 
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(5.13),  complex  periodicity  conditions  can  be  used  to  represent  the  effects  of  neighbor¬ 
ing  blade  passages  for  each  interblade  phase  (IBP)  angle  [121].  All  cascade  linearized 
CFD  computations  are  therefore  performed  in  the  frequency  domain  on  a  single  blade 
passage.  Similarly,  this  periodicity  can  be  exploited  to  provide  efficient  implementa¬ 
tions  for  creating  the  reduced-order  model. 

For  cascade  flows,  experimental  data  for  a  number  of  standard  configurations 
are  available  [122],  The  first  standard  configuration  with  a  steady-state  inflow  Mach 
number  of  0.18  and  flow  angle  (3  of  62°  is  considered.  The  cascade  operates  in  unsteady 
rigid  pitching  motion.  Both  experimental  data  and  other  CFD  data  [123, 124]  are 
taken  from  Boles  and  Fransson  [122]  and  Fransson  and  Verdon  [125].  Figures  5-2 
and  5-3  show  the  magnitude  and  the  phase  of  the  unsteady  pressure  coefficients  on 
the  first  blade  as  a  function  of  pitching  frequency  at  an  IBP  of  a  —  45°.  Figures 
5-4  and  5-5  show  the  pressure  coefficients  at  a  —  —45°.  It  can  be  seen  that  the  DG 
linearized  results  are  comparable  to  other  CFD  results  and  are  acceptably  close  to  the 
experimental  data.  It  should  be  noted  that  the  DG  method  is  very  sensitive  to  the 
geometry  representation  of  the  blade.  A  very  coarse  cascade  geometry  is  available 
from  Boles  and  Fransson  [122],  This  geometry  is  then  splined  to  obtain  smoother 
blade  surfaces.  It  is  expected  that  a  more  accurate  result  could  be  obtained  with  the 
DG  method  if  more  accurate  geometry  data  were  available. 


5.4  Geometric  Variability  Model 

Mistiming,  or  blade-to-blade  variation,  is  an  important  consideration  for  aeroelastic 
analysis  of  bladed  disks,  since  even  small  variations  among  blades  can  have  a  large 
impact  on  the  forced  response  and  consequently  the  high-cycle  fatigue  properties  of 
the  engine.  The  effects  of  blade  structural  mistiming  (variations  in  mass  and  stiffness 
properties)  have  been  extensively  studied,  see  for  example  Refs.  14-19;  however,  due 
to  the  prohibitively  high  computational  cost  of  performing  probabilistic  analysis  with 
a  CFD  model,  the  aerodynamic  effects  due  variations  in  geometry  are  less  under¬ 
stood.  Lim  et  al.  [126]  have  incorporated  geometric  mistiming  effects  into  structural 
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Figure  5-2:  First  standard  configuration  in  unsteady  pitching  motion  with  M  =  0.18, 
/ 3  =  62°.  Magnitude  and  phase  of  the  unsteady  pressure  coefficient  distribution  on 
the  lower  surface  with  a  =  45°. 


Figure  5-3:  First  standard  configuration  in  unsteady  pitching  motion  with  M  =  0.18, 
(3  =  62°.  Magnitude  and  phase  of  the  unsteady  pressure  coefficient  distribution  on 
the  upper  surface  with  er  =  45°. 
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Figure  5-4:  First  standard  configuration  in  unsteady  pitching  motion  with  M  =  0.18, 
/ 3  =  62°.  Magnitude  and  phase  of  the  unsteady  pressure  coefficient  distribution  on 
the  lower  surface  with  a  =  —45°. 


Figure  5-5:  First  standard  configuration  in  unsteady  pitching  motion  with  M  =  0.18, 
/ 3  =  62°.  Magnitude  and  phase  of  the  unsteady  pressure  coefficient  distribution  on 
the  upper  surface  with  a  =  —45°. 
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responses  of  bladed  disks.  In  their  work,  the  mode-acceleration  method  was  used  to 
convert  the  effect  of  geometric  mistiming  to  that  of  external  forces  of  the  tuned  disks. 
Truncated  sets  of  tuned  system  modes  compensated  by  static  modes — generated  by 
external  forces  that  were  constructed  from  mistiming — were  then  used  to  obtain  effi¬ 
cient  and  accurate  structural  reduced  models. 

Since  the  manufactured  geometric  mistiming  space  is  large,  Garzon  and  Darmofal 
[127-129],  Brown  et  al.  [130],  Ghiocel  [131],  and  Sinha  et  al.  [132]  have  used  the 
principle  component  analysis  (PGA)  [133]  to  construct  a  reduced  geometric  variability 
model.  It  was  found  that  a  handful  of  PGA  geometric  variability  modes  can  capture 
the  manufactured  variability  well.  In  particular,  Garzon  and  Darmofal  [127-129]  have 
used  the  reduced  geometric  variability  model  to  investigate  the  impact  of  geometric 
variability  on  axial  compressor  steady  aerodynamic  performance  using  Monte  Carlo 
simulation  based  on  a  large-scale  nonlinear  CFD  model.  They  found  that  the  mean 
loss  under  the  presence  of  geometric  mistiming  was  approximately  20%  larger  than 
the  nominal  loss.  Since  each  large-scale  nonlinear  CFD  simulation  is  very  expensive, 
parallel  computers  were  used  in  order  to  perform  the  probabilistic  analysis  using 
Monte  Carlo  simulation.  Here  we  consider  incorporating  the  effects  of  geometry 
variability  into  the  linearized  unsteady  CFD  model. 

Following  Ref.  127,  a  general  geometry,  g,  can  be  expressed  as 

n3 

9  =  9n  +  g  +  ^  (5.14) 

i— 1 

where  gn  is  the  nominal  geometry,  g  is  the  average  geometric  variation,  vt  are  ge¬ 
ometric  mode  shapes,  and  ns  is  the  number  of  mode  shapes  used  to  represent  the 
variation  in  geometry.  The  geometric  mode  shapes  could  be  computed,  for  example, 
by  performing  the  PCA  on  a  manufacturing  sample  of  geometries.  In  that  case,  the 
parameters  zt  in  (5.14)  are  random  numbers  normally  distributed  with  zero  mean 
and  unity  variance,  Zi  G  IV(0,  l)jj  and  a,  is  the  standard  deviation  of  the  geometric 
data  attributable  to  the  ith  mode;  thus  the  product  cr^  is  the  amount  by  which  the 
mode  Vi  contributes  to  the  geometry  g.  A  detailed  description  of  the  methodology 
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underlying  this  geometric  model  can  be  found  in  Ref.  127. 

The  key  assumption  in  (5.14)  is  that  any  manufacturing  geometry  is  a  linear 
combination  of  geometries  in  the  manufacturing  geometry  sample.  In  other  words, 
the  manufacturing  geometry  sample  is  assumed  to  be  large  enough  to  span  the  entire 
mistiming  space. 

5.5  Linearized  Unsteady  CFD  Model  with  Geo¬ 
metric  Variability 

Using  the  model  (5.14),  a  general  geometry  g( z)  is  specified  by  the  parameter  vector 
z  =  [z\,  Z2,  ■  •  • ,  zns]T,  which  describes  the  geometry  variability  in  terms  of  the  ge¬ 
ometry  modes.  The  linearized  CFD  system  corresponding  to  geometry  g( z)  is  given 
by 


E(wm(0(z)),0(z))x  =  A(wm(0(z)),0(z))x  +  B(wm(0(z)),0(z))u,  (5.15) 

y  =  C(wss(<?(z)),  g(z))x,  (5.16) 

where  the  CFD  system  matrices  E.  A,  B  and  C  are  in  general  both  a  function  of 
the  geometry,  g( z),  and  of  the  steady-state  solution,  wss(g(z)),  which  is  itself  also  a 
function  of  the  geometry.  To  solve  the  CFD  system  (5.15),  (5.16),  for  each  geometry 
g  we  must  firstly  compute  the  steady-state  solution,  wss(g(z)),  secondly  evaluate  the 
linearized  matrices  E,  A,  B  and  C,  and  thirdly  solve  the  resulting  large-scale  linear 
system.  This  is  a  computationally  prohibitive  proposition  for  applications  such  as 
probabilistic  analysis,  where  thousands  of  geometry  perturbations  may  be  analyzed 
over  many  random  samples  z.  For  example,  if  one  such  analysis  takes  three  minutes 
to  perform,  then  50,000  analyses  would  take  more  than  three  months  of  CPU  time! 

For  convenience  of  notation,  we  write  the  dependence  of  the  CFD  matrices  on 
the  parameter  z  as  E(wss(g(z)),  g(z))  =  E(z),  A(wss(#(z)),  g(z))  =  A(z), 

B(wss(g(z)),  g(z))  =  B(z),  and  C(wss(g(z)),  g(z))  =  C(z).  We  use  the  expan¬ 

sion  given  by  equation  (5.14),  which  represents  a  general  geometry  as  a  perturbation 
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about  the  average  geometry  g0  =  gn  +  g,  to  derive  an  approximate  model  for  repre¬ 
senting  the  effects  of  geometry  variations.  Instead  of  computing  the  linearized  CFD 
matrices  exactly  for  any  random  variability  z,  we  choose  to  linearize  the  relationships 
E(z),  A(z),  B(z),  and  C(z)  [65,66].  A  more  general  approach  for  a  general  nonlinear 
function  can  be  found  in  Barrault  et  al.  [61,62,68].  We  define  the  linearized  unsteady 
CFD  model  for  the  average  geometry  g0  =  gn  +  g  by  the  matrices  E0,  A0,  B0,  and 
Co,  with  corresponding  solution  x0.  That  is,  for  z  =  0  we  have 


EqXo  —  A0x0  +  B0u, 


(5.17) 


y0  =  c0x0. 


(5.18) 


Using  a  Taylor  series  expansion  about  z  =  0  for  the  matrix  A(z)  gives 


f)  A 

A(z,^A„+- 


Z\  +  . . .  + 


dA 


z=0 


dZn 


+  *  *  *  > 


z=0 


(5.19) 


where  the  matrix  partial  derivatives  denote  componentwise  derivatives,  which  can  be 
evaluated  through  application  of  the  chain  rule.  These  derivatives  are  evaluated  at 
average  conditions,  z  =  0.  The  matrices  E(z),B(z)  and  C(z)  can  be  expanded  using 
formulae  analogous  to  (5.19). 


If  the  geometric  variability  (given  by  the  product  cr^)  is  sufficiently  small,  the 
constant  and  linear  terms  in  the  Taylor  expansion  (5.19)  are  sufficient  to  approximate 
the  linearized  matrices  A(z)  accurately,  that  is, 


A(z)  «  An 


dA 

dz\ 


Z\  +  . . .  + 


dA 


z=0 


dzn 


z=0 


(5.20) 


For  i  —  1,2, ,  ns,  we  define 
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dE 

dzi 


A, 


z=0 


dA 
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Z  — 0 


dB 
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Ci 


z=0 


d  C 


dzi 


1 

z=0 


(5.21) 


where  the  matrices  E,,  A,,  Bj,  and  C.,  can  be  computed,  for  example,  using  a  finite 
difference  approximation  of  the  respective  derivatives.  The  approximate  linearized 


118 


CFD  model  for  any  geometric  variability  z  is  then  given  by 


(5.22) 

(5.23) 


It  should  be  noted  here  that  a  number  of  large-scale  steady  state  CFD  solves  are 
required  in  order  to  determine  the  matrices  A0,  B0,  Co,  E0,  A,:,  Bj,  C,  and  Ej.  For 
example,  if  central  difference  approximations  to  the  matrix  derivatives  are  used,  a 
total  of  2  Y^ili  +1  large-scale  steady  state  CFD  solves  is  required.  This  is  a  one¬ 
time  offline  cost;  once  the  matrices  are  computed,  the  approximate  linearized  system 
(5.22),  (5.23)  can  be  readily  evaluated  for  an  arbitrary  geometry  g( z)  without  running 
the  CFD  steady  solver.  It  is  important  to  know  that  the  size  and  components  of  A(z) 
are  a  function  of  the  CFD  grid.  The  CFD  grid  is  in  turn  a  function  of  the  geometry, 
which  is  in  turn  a  function  of  geometric  variability  z.  This  implies  that,  for  the  finite 
difference  to  be  accurate,  the  grid  generation  must  satisfy  the  requirements  that  the 
size  of  the  linearized  matrices  be  the  same  for  any  geometric  variability  and  that 
the  components  of  the  linearized  matrices  be  a  smooth  function  of  the  geometric 
variability  z.  In  order  to  satisfy  these  requirements,  we  first  generate  a  CFD  grid 
for  the  average  geometry  g0  =  gn  +  g.  The  grid  for  any  new  geometry  corresponding 
to  a  nonzero  geometric  variability  z  is  then  generated  by  adding  the  perturbation 
aizivt  to  the  boundary  grid  points  and  computing  the  new  grid  points  using  the 
Jacobi  smoothing  in  (5.10). 


It  should  also  be  noted  that  the  model  (5.22),  (5.23)  is  valid  only  for  small  varia¬ 
tions  from  the  average  geometry.  Larger  variations  will  incur  larger  errors,  clue  to  the 
neglect  of  the  higher-order  terms  in  the  Taylor  series  expansion.  In  the  next  section, 
we  present  some  analyses  to  quantify  these  errors.  Even  with  this  restriction,  the 
model  is  useful  for  many  applications  where  small  geometric  variations  are  of  inter- 
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est;  however,  the  approximate  linearized  model  is  still  of  high  dimension,  and  thus  is 
computationally  too  expensive  for  applications  such  as  probabilistic  analysis  in  which 
one  needs  to  determine  the  unsteady  aerodynamic  response  for  many  random  geome¬ 
tries.  In  Chapter  6  we  further  reduce  the  cost  of  solving  the  approximate  linearized 
system  by  using  the  model  reduction  method  developed  in  Chapter  3  to  create  a 
reduced-order  model  that  is  accurate  over  both  time  and  the  geometric  parameter 
space,  described  here  by  the  vector  z. 


5.6  Linearized  Geometric  Variability  CFD  Model 
Validation 

Results  are  presented  for  forced  response  of  a  subsonic  rotor  blade  that  moves  in 
unsteady  rigid  motion.  The  flow  is  modeled  using  the  two-dimensional  Euler  equations 
written  at  the  blade  mid-section.  The  average  geometry  of  the  blade  is  shown  in 
Figure  5-6  along  with  the  unstructured  grid  for  a  single  blade  passage,  which  contains 
4292  triangular  elements.  The  Euler  equations  are  discretized  in  space  with  the 
discontinuous  Galerkin  (DG)  method  described  in  Section  5.1.  For  the  case  considered 
here,  the  incoming  steady-state  flow  has  a  Mach  number  of  M  —  0.113  and  a  flow  angle 
of  (3  —  59°.  Flow  tangency  boundary  conditions  are  applied  on  the  blade  surfaces. 
Since  the  rotor  is  cyclically  symmetric,  the  steady  flow  in  each  blade  passage  is  the 
same  and  the  steady-state  solution  can  be  computed  on  a  computational  domain  that 
describes  just  a  single  blade  passage.  Periodic  boundary  conditions  are  applied  on  the 
upper  and  lower  boundaries  of  the  grid  to  represent  the  effects  of  neighboring  blade 
passages. 

A  linearized  model  is  derived  for  unsteady  flow  computations  by  assuming  that 
the  unsteady  flow  is  a  small  deviation  from  steady  state  as  described  in  Section  5.2. 
An  affine  dependence  of  the  linearized  system  matrices  on  the  blade  geometries  is 
derived  using  the  method  described  in  Section  5.5.  This  leads  to  a  system  of  the 
form  (5.22),  (5.23),  where  the  state  vector,  x(f),  contains  the  unknown  perturbation 
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Figure  5-6:  Geometry  and  CFD  mesh  for  a  single  blade  passage. 


flow  quantities  (density,  Cartesian  momentum  components  and  energy).  For  the  DG 
formulation,  the  states  are  the  coefficients  corresponding  to  each  nodal  finite  element 
shape  function.  Using  linear  elements,  there  are  12  degrees  of  freedom  per  element, 
giving  a  total  state-space  size  of  n  =  51,  504  states  per  blade  passage.  For  the  problem 
considered  here,  the  forcing  input,  u(f),  describes  the  unsteady  motion  of  each  blade, 
which  in  this  case  is  assumed  to  be  rigid  plunging  motion  (vertical  motion  with  no 
rotation).  The  outputs  of  interest,  y (f),  are  the  unsteady  lift  forces  and  pitching 
moments  generated  on  each  blade.  The  initial  perturbation  flow  is  given  by  x°  =  0. 

Geometric  modes  were  computed  using  a  PCA  model  of  data  from  145  actual 
blades,  measured  at  thirteen  sections  along  the  radial  direction.  The  mid-section 
geometries  were  then  extracted.  Thus  the  parameter  vector  z  contains  the  normally 
distributed  random  variables  that  describe  perturbations  in  the  geometry  of  each 
blade  according  to  the  model  (5.14).  Since  the  approximate  linearized  CFD  model  is 
only  valid  for  small  variations  from  the  average  geometry,  the  standard  deviation  of 
the  actual  manufacturing  data  was  reduced  by  a  factor  of  6.  As  the  results  below  show, 
this  ensures  that  the  geometric  model  remains  in  its  region  of  applicability;  however, 
it  also  highlights  a  limitation  in  the  geometric  model  used  here.  By  including  ideas 
from  [61,62,68]  to  handle  a  general  nonlinear  term  in  an  efficient  way  together  with 
the  framework  proposed  here,  a  more  general  geometric  model  could  be  derived  that 
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is  applicable  for  larger  geometric  deviations. 

In  Figure  5-7,  we  consider  a  geometric  model  that  uses  the  two  dominant  variabil¬ 
ity  modes,  ns  —  2.  The  hgure  shows  the  lift  coefficient,  Cl,  and  moment  coefficient, 
Cm,  of  a  blade  in  response  to  a  pulse  input  in  plunge  for  a  particular  geometry 
that  corresponds  to  z\  =  1.59,  z2  =  1.59.  The  response  is  computed  using  the  exact 
linearized  CFD  model,  i.e.  the  system  (5.15),  (5.16)  and  the  approximate  linearized 
model  (5.22),  (5.23)  with  ns  =  2  geometry  modes.  For  reference,  the  response  of 
the  nominal  blade  is  also  shown  in  the  hgure.  It  can  be  seen  that  despite  the  small 
perturbation  in  geometry,  the  change  in  lift  and  moment  coefficient  responses  is  sig¬ 
nificant.  The  approximate  linearized  geometric  model  captures  the  unsteady  response 
accurately. 


Figure  5-7:  Lift  coefficient,  Cl,  and  moment  coefficient,  Cm,  in  response  to  a  pulse 
input  in  blade  plunge  displacement  for  the  nominal  geometry  and  a  perturbed  geom¬ 
etry  described  by  two  geometric  PCA  modes  with  Z\  =  1.59,  z2  =  1.59.  Perturbed 
geometry  results  are  computed  with  both  the  exact  and  approximate  linearized  CFD 
model. 


Table  5.1  shows  the  error  in  lift  and  moment  outputs  due  to  the  linearized  geome¬ 
try  approximation  for  several  different  blade  geometries  with  a  pulse  input  in  plunge. 
The  error  e  is  defined  as  the  2-norm  of  the  difference  between  the  approximate  and 
the  exact  linearized  output  as  a  percentage  of  the  change  between  the  exact  and  the 
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nominal  output, 


\[$  \\Ye  ~  YaWldt 

sif  wye-yoWidt 

where  ye,  ya,  and  y0  are  respectively  the  exact,  approximate,  and  nominal  outputs. 
In  the  table,  ecM  denotes  the  error  in  moment  coefficient  response,  while  ecL  denotes 
the  error  in  lift  coefficient  response.  In  general,  we  expect  the  quality  of  the  approx¬ 
imate  model  to  be  compromised  as  the  size  of  the  geometric  perturbation  increases. 
The  errors  shown  in  Table  5.1  for  blade  geometries  in  the  tails  of  the  distribution, 
i.e.  those  with  large  geometry  variation,  are  deemed  to  be  acceptable  for  the  proba¬ 
bilistic  application  of  interest  here  (although  again  we  note  that  the  actual  variations 
observed  in  manufacturing  data  are  larger).  In  fact,  as  will  be  shown  in  probabilistic 
analysis  in  Chapter  6,  the  error  in  aggregate  quantities  will  be  shown  to  be  less.  For 
applications  where  greater  accuracy  for  large  geometry  variations  is  important  (for 
example,  determining  the  probability  of  failure  would  require  the  tail  of  the  distribu¬ 
tion  to  be  resolved  accurately),  the  results  suggest  that  the  approximate  linearized 
CFD  system  is  not  appropriate.  In  such  cases,  one  might  consider  including  high 
order  terms  in  (5.20),  the  Taylor  series  expansion  of  the  CFD  matrices. 


x  100%, 


(5.24) 


Table  5.1:  Error  in  approximate  linearized  model  predictions  for  a  pulse  input  in 
blade  displacement  for  several  different  geometries. 


Variability  amplitudes 

ecM(%) 

ecL  (%) 

Z!  =  1.59,  *2  =  1.59 

5.04 

2.6 

=  1.59,  z2  =  -1.59 

0.3 

0.1 

Z!  =  —1.59,  z2  =  -1.59 

2.0 

0.8 

Zi  =  3.0,  z2  =  3.0 

16.6 

9.2 

z\  =  3.0,  z2  =  —3.0 

4.1 

2.3 

z\  =  —3.0,  z2  =  —3.0 

12.4 

4.7 
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Chapter  6 


CFD  Probabilistic  Analysis 
Application 


The  chapter  begins  with  a  numerical  result  for  the  Petrov- Galerkin  projection  in  Sec¬ 
tion  2.3.  The  Petrov-Galerkin  projection  is  then  applied  to  obtain  reduced  models 
for  all  examples  that  follow.  The  model-constrained  greedy-based  adaptive  sampling 
approach  developed  in  Chapter  3  is  applied  to  probabilistic  analysis  problems  that  de¬ 
pend  on  a  large  number  of  parameters.  In  particular,  probabilistic  analysis  problems 
with  four  and  ten  geometric  variability  parameters  will  be  considered.  Finally,  we 
compare  the  model-constrained  sampling  method  with  statistically-based  sampling 
methods  in  the  context  of  snapshot  generation  for  model  reduction  for  the  probabilis¬ 
tic  analysis  problems  with  four  and  ten  geometric  variability  parameters. 


6.1  Galerkin  Projection  versus  Petrov-Galerkin  Pro¬ 
jection 

First,  we  present  results  to  emphasize  the  importance  of  using  an  appropriate  pro¬ 
jection  basis  to  perform  the  reduction.  In  particular,  we  verify  the  stability  and  the 
convergence  of  the  projection  method  discussed  in  Section  2.3,  in  which  the  resid¬ 
uals  are  minimized  sequentially.  We  consider  a  subsonic  rotor  blade  that  moves  in 
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unsteady  rigid  motion,  as  described  in  Section  5.6.  In  this  example,  the  nominal 
geometry  with  zero  interblade  phase  angle,  i.e.  IBP  =  0°,  is  used.  The  linearized 
CFD  system  (5.11)  therefore  does  not  depend  on  the  geometric  variability  vector  z, 
and  the  unsteady  flow  can  be  computed  on  a  computational  domain  which  describes 
just  a  single  blade  passage.  The  forcing  input,  u(£)  is  assumed  to  be  rigid  plunging 
motion  (vertical  motion  with  no  rotation).  The  output  of  interest,  y,  is  the  unsteady 
lift  force  generated  on  the  blade.  The  initial  perturbation  flow  is  given  by  x°  =  0. 


Snapshots  are  taken  by  computing  the  response  of  the  blade  to  a  pulse  input  in 
plunging  motion.  For  this  input,  the  blade  vertical  position  as  a  function  of  time  is 
given  by 

h(t)  =  he~9{t-to)2,  (6.1) 

where  the  parameters  h  —  0.1,  g  —  0.02,  and  t0  =  40  are  chosen  based  on  the  range 
of  motions  that  are  expected  in  practice,  and  all  quantities  are  non-dimensionalized 
with  the  blade  chord  as  a  reference  length  and  the  inlet  speed  of  sound  as  a  reference 
velocity.  The  unsteady  simulation  is  performed  with  a  timestep  of  At  =  0.1  from 
t  —  0  to  tf  —  200.  A  set  of  POD  basis  vectors  is  computed  from  this  collection  of 
2000  snapshots  and  is  used  as  the  trial  reduced  basis  <h. 


This  example  is  of  interest  since  the  reduced  model  is  unstable,  depending  on 
the  number  of  reduced  basis  vectors,  if  the  usual  Galerkin  projection  is  used.  Here 
we  guarantee  the  stability  of  the  reduced  model  by  minimizing  the  residual  at  each 
time  step  sequentially,  and  hence  by  the  Petrov- Galerkin  projection  (2.35)  in  Chap¬ 
ter  2.  The  comparison  between  the  Galerkin-projection-based  and  Petrov-Galerkin- 
projection-based  reduced  models  can  be  seen  in  Figure  6-1.  Note  that  we  have  used 
the  same  POD  basis  vectors  for  both  methods.  As  expected,  the  Petrov- Galerkin  ap¬ 
proach  is  stable  and  improves  the  reduced  model  as  more  basis  vectors  are  used  while 
the  Galerkin  approach  yields  unstable  models  for  any  basis  size  less  than  m  =  10. 
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time  time 


(a)  Galerkin  POD  with  2  basis  vectors  (b)  Petrov-Galerkin  POD  with  2  basis  vec¬ 
tors 


time  time 


(c)  Galerkin  POD  with  8  basis  vectors  (d)  Petrov-Galerkin  POD  with  8  basis  vec¬ 
tors 


Figure  6-1:  Comparison  between  the  conventional  Galerkin  projection  and  the  Petrov- 
Galerkin  projection  approaches.  Solid  lines  are  the  full  (exact)  solution  and  dashed 
lines  are  reduced  model  results. 
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6.2  Probabilistic  Analysis  Application 


The  model-constrained  adaptive  model  reduction  method  is  applied  to  probabilistic 
analysis  of  a  subsonic  rotor  blade  that  moves  in  unsteady  rigid  motion.  Using  Monte 
Carlo  simulation  (MCS)  of  a  CFD  model  to  quantify  the  impact  of  geometric  vari¬ 
ability  on  unsteady  performance  is  a  computationally  prohibitive  proposition.  For 
example,  if  the  unsteady  analysis  for  one  geometry  takes  three  minutes  to  compute 
(a  conservative  estimate),  the  0(50,000)  such  analyses  that  would  be  required  for  a 
MCS  would  take  more  than  three  months  of  CPU  time.  Therefore,  we  desire  to  ob¬ 
tain  a  reduced-order  model  that  captures  both  unsteady  response  and  variation  over 
blade  geometries.  Our  method  combines  the  reduced  geometric  variability  model 
and  the  model-constrained  adaptive  sampling  procedure  of  Algorithm  3.1  to  obtain 
a  reduced-order  model  that  is  valid  over  a  range  of  forcing  frequencies,  aerodynamic 
damping,  and  small  perturbations  in  blade  geometries,  and  thus  enables  fast  and 
accurate  probabilistic  analysis. 

Results  are  shown  here  for  the  case  of  two  blades  moving  with  an  interblade  phase 
angle  of  180°.  First,  each  blade  geometry  is  represented  by  two  variability  modes, 
giving  d  —  ns  —  4  geometric  parameters  in  this  example.  Applying  the  adaptive 
model  reduction  methodology  with  e  =  10-4  and  with  the  lift  coefficients  of  the 
blades  as  the  outputs  of  interest  yields  a  reduced-order  model  of  size  nr  =  201  (for 
two  blades).  Again,  in  each  greedy  cycle,  the  number  of  Newton  steps  is  observed 
to  scale  as  O(d).  Algorithm  3.1  requires  21  greedy  cycles,  over  which  a  total  of  21 
optimal  parameter  points  are  found.  Recall  that  when  using  the  error  indicator,  the 
snapshots  are  only  computed  at  the  optimal  parameters.  Therefore,  the  full  model 
is  solved  21  times,  and  this  computational  cost  dominates  the  other  calculations.  At 
each  optimal  parameter  point,  unlike  the  steady  cases  in  which  only  one  snapshot  is 
found,  for  unsteady  cases,  401  snapshots  are  taken  uniformly  over  the  time  horizon 
tf  —  200.  These  snapshots  are  then  used  to  update  the  reduced  basis  using  a  strategy 
similar  to  the  third  method  in  Table  4.3.  That  is,  the  newly  computed  snapshots  are 
pre-processed  by  the  POD  method  (where  rj  =  0.99999999),  and  the  dominant  POD 
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basis  vectors  are  then  added  to  the  current  reduced  basis  using  the  Gram-Schmidt 
procedure.  The  geometric  variability  parameters,  z,  are  random  numbers  normally 
distributed  with  zero  mean  and  unity  variance,  99.7%  of  which  are  distributed  in 
the  interval  [—3.0, 3.0].  The  bound  constraints  for  all  parameters,  in  the  greedy 
optimization  problem,  are  therefore  chosen  to  be 

—3.0  <  Zi  <  3.0,  %  —  1, _ ,  d.  (6.2) 

Finally,  the  training  input  u(f)  is  chosen  to  be  the  pulse  function  in  (6.1). 

We  now  have  a  reduced  model  of  size  nr  =  201  that  accurately  captures  the 
unsteady  response  of  the  original  two-blade  system  with  n  =  103,  008  states  over  the 
range  of  geometries  described  by  the  four  geometric  parameters.  As  an  example  of  an 
application  for  which  this  reduced  model  is  useful,  we  consider  probabilistic  analysis 
of  the  system.  Specifically,  we  consider  the  impact  of  blade  geometry  variabilities  on 
the  work  per  cycle,  which  is  defined  as  the  integral  of  the  blade  motion  times  the  lift 
force  over  one  unsteady  cycle.  A  MGS  is  performed  in  which  10,000  blade  geometries 
are  selected  randomly  from  the  given  distributions  for  each  blade.  The  same  10,000 
geometries  are  analyzed  using  the  approximate  linearized  CFD  model  and  the  reduced 
model.  Figure  6-2  shows  the  resulting  probability  density  functions  (PDFs)  of  work 
per  cycle  for  the  first  blade,  computed  using  the  approximate  linearized  CFD  model 
and  the  reduced-order  model.  Figure  6-3  shows  the  PDFs  of  work  per  cycle  for  the 
second  blade.  Table  6.1  shows  that  the  online  CPU  time  required  to  compute  the 
reduced  model  MGS  is  a  factor  of  2414  times  smaller  than  that  required  for  the  CFD 
MGS.  As  can  be  also  seen,  the  savings  in  the  online  cost  are  substantial,  and  more 
than  justify  the  offline  cost  required  to  compute  the  reduced  model.  In  practice,  many 
more  than  10,000  blade  geometries  are  required  to  obtain  a  converged  MGS;  in  this 
case,  the  computational  cost  of  using  the  CFD  model  becomes  prohibitive.  These 
computational  results  were  obtained  on  a  dual  core  64-bit  personal  computer  with 
3.2GHz  Pentium  processor.  Tables  6.1  also  shows  the  number  of  nonzeros  in  the  full 
system  matrix  (sparse),  A,  and  the  reduced  system  matrix  (dense),  Ar. 
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(a)  Full  WPC  for  Blade  1  (b)  Reduced  WPC  for  Blade  1 

Figure  6-2:  Comparison  of  linearized  CFD  and  reduced-order  model  predictions  of 
work  per  cycle  (WPC)  for  Blade  1.  MCS  results  are  shown  for  10,000  blade  geometries 
with  four  parameters.  The  same  geometries  are  analyzed  in  each  case.  Dashed  line 
denotes  the  mean. 
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(a)  Full  WPC  for  Blade  2 


(b)  Reduced  WPC  for  Blade  2 


Figure  6-3:  Comparison  of  linearized  CFD  and  reduced-order  model  predictions  of 
work  per  cycle  for  Blade  2.  MCS  results  are  shown  for  10,000  blade  geometries  for  the 
case  of  four  parameters.  The  same  geometries  were  analyzed  in  each  case.  Dashed 
line  denotes  the  mean. 
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Table  6.1  also  compares  the  statistics  of  the  two  distributions.  It  can  be  seen  from 
Figure  6-2,  Figure  6-3  and  Table  6.1  that  the  reduced-order  model  predicts  the  mean, 
variance  and  shape  of  the  distribution  of  work  per  cycle  accurately.  To  further  verify 
the  quality  of  the  reduced  model,  we  apply  the  Kolmogorov- Smirnov  method  [134], 
to  test  whether  the  reduced  work  per  cycle  results  and  the  full  work  per  cycle  results 
are  drawn  from  a  same  distribution.  The  results  show  that  we  cannot  reject  the 
hypothesis  that  the  distribution  is  the  same  at  a  5%  significance  level.  It  should  be 
pointed  out  the  Kolmogorov-Smirnov  method  is  based  on  the  maximum  difference 
between  the  empirical  cumulative  distribution  functions  of  work  per  cycle  of  the  full 
and  the  reduced  models.  As  shown  in  Figure  6-4  for  Blade  1,  the  empirical  cumulative 
distribution  functions  of  work  per  cycle  of  the  full  and  the  reduced  models  are  close 
to  one  another. 

To  further  compare  the  reduced-order  and  CFD  results,  we  pick  four  particular 
geometries  corresponding  to  the  left  tail,  right  tail,  mid-left  and  mid-right  locations 
on  the  PDF  of  the  first  blade  as  indicated  by  the  circles  in  Figure  6-2(a).  In  Table  6.2, 
the  work  per  cycle  is  given  for  these  four  blade  geometries  as  computed  by  the  exact 
CFD  model,  the  approximate  linearized  CFD  model,  and  the  reduced-order  model. 
The  table  shows  that  again  the  approximate  linearized  CFD  is  in  good  agreement 
with  the  exact  CFD,  especially  for  the  mid-left  and  mid-right  cases,  which  have 


Table  6.1:  Linearized  CFD  model  and  reduced-order  model  MCS  results  for  the  case 
of  four  parameters.  Work  per  cycle  (WPG)  is  predicted  for  blade  plunging  motion  at 
an  interblade  phase  angle  of  180°  for  10,000  randomly  selected  blade  geometries. 


CFD 

Reduced 

Model  size 

103,008 

201 

Number  of  nonzeros 

2,846,056 

40,401 

Offline  cost 

— 

2.8  hours 

Online  cost 

501.1  hours 

0.21  hours 

Blade  1  WPC  mean 

-1.8572 

-1.8573 

Blade  1  WPC  variance 

2.687e-4 

2.6819e-4 

Blade  2  WPC  mean 

-1.8581 

-1.8580 

Blade  2  WPC  variance 

2.797e-4 

2.799e-4 
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Empirical  CDF 


Figure  6-4:  Comparison  of  the  approximate  linearized  CFD  and  the  reduced-order 
model  empirical  cumulative  distribution  functions  (CDF)  of  work  per  cycle  (WPC). 
The  results  are  shown  for  10,000  blade  geometries  for  the  case  of  four  parameters. 
The  same  geometries  were  analyzed  for  both  the  full  CFD  and  its  reduced  model. 


smaller  variability.  In  addition,  the  effectiveness  of  the  adaptive  model  reduction 
methodology  of  Algorithm  3.1  can  be  seen  from  the  good  agreement  between  the 
approximate  linearized  CFD  and  the  reduced  results. 


Table  6.2:  Exact  CFD,  approximate  CFD,  and  reduced-order  model  work  per  cycle 
prediction  for  the  four  geometries  indicated  in  Figure  6-2 (a). 


Exact 

Approximate 

Reduced 

Left  tail 

-1.8973 

-1.9056 

-1.9060 

Mid-left 

-1.8637 

-1.8636 

-1.8638 

Mid-right 

-1.8459 

-1.8455 

-1.8458 

Right  tail 

-1.8014 

-1.8086 

-1.8088 

Next,  we  consider  the  case  in  which  each  blade  geometry  is  represented  by  five 
variability  modes,  giving  d  =  ns  =  10  geometric  parameters  in  this  example.  Applying 
the  adaptive  model  reduction  methodology  with  e  =  10-4  and  with  the  lift  coefficients 
of  the  blades  as  the  outputs  of  interest  yields  a  reduced-order  model  of  size  nr  =  290 
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(for  two  blades).  Algorithm  3.1  requires  29  greedy  cycles,  over  which  a  total  of  29 
optimal  parameter  points  are  found. 


We  now  have  a  reduced  model  of  size  nr  =  290  that  accurately  captures  the 
unsteady  response  of  the  original  two-blade  system  with  n  =  103,  008  states  over 
the  range  of  geometries  described  by  the  ten  geometric  parameters.  A  MCS  is  then 
performed  in  which  10,000  blade  geometries  are  selected  randomly  from  the  given 
distributions  for  each  blade.  The  same  10,000  geometries  are  analyzed  using  the 
approximate  linearized  CFD  model  and  the  reduced  model.  Figure  6-5  shows  the 
PDFs  of  work  per  cycle  for  the  first  blade.  Figure  6-6  shows  the  PDFs  of  work  per 
cycle  for  the  second  blade.  Table  6.3  shows  that  the  online  CPU  time  required  to 
compute  the  reduced  model  MCS  is  a  factor  of  468  times  smaller  than  that  required 
for  the  CFD  MCS.  Note  that  the  observed  speed-up  factor  in  this  case  is  smaller  than 
that  observed  in  the  case  of  four  parameters.  This  is  due  to  the  fact  that  not  only 
we  have  more  parameters  but  also  the  reduced  model  size  is  now  larger.  Nonetheless, 
the  savings  in  the  online  cost  are  substantial  and  offset  the  offline  cost. 


JL 

-2.5  -2  -1.5  -1 

WPC 

(a)  Full  WPC  for  Blade  1 

Figure  6-5:  Comparison  of  linearized  CFD  and  reduced-order  model  predictions  of 
work  per  cycle  for  Blade  1.  MCS  results  are  shown  for  10,000  blade  geometries  with 
ten  parameters.  The  same  geometries  are  analyzed  in  each  case.  Dashed  line  denotes 
the  mean. 

Table  6.3  also  compares  the  statistics  of  the  two  distributions.  It  can  be  seen 
from  Figure  6-5,  Figure  6-6  and  Table  6.3  that  the  reduced-order  model  predicts  the 
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(a)  Full  WPC  for  Blade  2  (b)  Reduced  WPC  for  Blade  2 


Figure  6-6:  Comparison  of  linearized  CFD  and  reduced-order  model  predictions  of 
work  per  cycle  for  Blade  2.  MCS  results  are  shown  for  10,000  blade  geometries  for 
the  case  of  ten  parameters.  The  same  geometries  are  analyzed  in  each  case.  Dashed 
line  denotes  the  mean. 


Table  6.3:  Linearized  CFD  model  and  reduced-order  model  MCS  results  for  the  case 
of  ten  parameters.  Work  per  cycle  (WPC)  is  predicted  for  blade  plunging  motion  at 
an  interblade  phase  angle  of  180°  for  10,000  randomly  selected  blade  geometries. 


CFD 

Reduced 

Model  size 

103,008 

290 

Number  of  nonzeros 

2,846,056 

84,100 

Offline  cost 

— 

10.92  hours 

Online  cost 

515.61  hours 

1.10  hours 

Blade  1  WPC  mean 

-1.8583 

-1.8515 

Blade  1  WPC  variance 

0.0503 

0.0506 

Blade  2  WPC  mean 

-1.8599 

-1.8583 

Blade  2  WPC  variance 

0.0136 

0.0138 
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mean,  variance  and  shape  of  the  distribution  of  work  per  cycle  accurately.  Again,  to 
further  verify  the  quality  of  the  reduced  model,  we  apply  the  Kolmogorov-Smirnov 
method,  to  test  whether  the  reduced  work  per  cycle  results  and  the  full  work  per 
cycle  results  are  drawn  from  a  same  distribution.  The  results  show  that  we  cannot 
reject  the  hypothesis  that  the  distribution  is  the  same  at  a  5%  significance  level. 

It  should  be  also  pointed  out  that  compared  to  the  case  with  four  parameters,  the 
variances  of  work  per  cycle  of  the  10-parameter  case  arc  larger;  however  the  means  are 
almost  the  same.  This  is  because  as  more  geometric  variability  parameters  are  added, 
and  hence  more  POD  modes  (which  are  variations  about  the  average  geometry)  are 
added  in  (5.14),  more  variability  is  added  to  the  blade  geometries.  It  is  therefore 
expected  that  the  variances  of  work  per  cycle  increase  as  more  geometric  variability 
parameters  are  included,  while  the  means  stay  the  same. 

Next,  we  compare  the  model-constrained  sampling  method  with  statistically-based 
sampling  methods  in  the  context  of  snapshot  generation  for  model  reduction  of  the 
linearized  CFD  model.  In  particular,  we  compare  our  model-constrained  sampling 
with  LHS  sampling,  logarithmic  random  sampling,  and  CVT  sampling.  Tables  6.4 
and  6.5  show  the  percentage  errors  in  predicting  the  means  and  variances  of  work 
per  cycle  of  the  resulting  reduced  models  for  the  case  with  four  and  ten  geometric 
variability  parameters.  The  percentage  error  in  predicting  the  quantity  h  ( h  can  be 
the  mean  or  the  variance  of  work  per  cycle)  is  defined  as 

Percentage  error  =  ^lcdu^cd — ^ful1^  x  100%,  (6.3) 

I  'huii  I 

where  hfuu  and  hreduced  are  the  corresponding  values  obtained  from  the  full  and  the 
reduced  models,  respectively.  As  can  be  seen,  the  quality  of  the  resulting  reduced 
models  is  similar.  While  the  system  matrices  of  the  steady  thermal  fin  problem  in 
(4.14)~(4.16)  are  nonlinear  functions  of  the  parameters,  the  system  matrices  of  the 
linearized  CFD  model  in  (5.22)-(5.23)  are  linear  functions  of  the  geometric  vari¬ 
ability  parameters.  As  a  result,  the  parameter-to-output  relation  for  the  linearized 
CFD  model  is  expected  to  be  weakly  nonlinear,  and  thus  statistically-based  sampling 
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methods  may  be  expected  to  perform  adequately.  In  particular,  for  problems  that 
are  linear  in  the  state  vector,  if  the  parameter-to-output  relation  is  also  linear,  with 
one  sample  point,  all  methods  yield  a  same  reduced  model.  In  addition,  as  mentioned 
in  Section  6.2,  the  bound  constraints  of  the  parameters  in  the  linearized  CFD  model 
are  relatively  small  as  compared  to  those  in  the  steady  thermal  hn  problem.  For 
a  small  parameter  range,  unless  the  problem  under  consideration  is  very  nonlinear 
with  respect  to  the  parameters  in  that  small  parameter  range,  it  is  expected  that  the 
sampling  methods  will  perform  adequately. 


Table  6.4:  Model-constrained  (MC)  adaptive  sampling  method  versus  logarithmic 
random  (LogRandom)  sampling,  LHS  and  CVT  for  the  case  with  four  geometric 
variability  parameters.  For  all  methods,  21  sample  points  are  generated,  and  the 
percentage  errors  in  predictions  of  means  and  variances  of  work  per  cycle  (WPC)  of 
the  resulting  reduced  models  are  shown  for  10,000  randomly  selected  blade  geometries. 


MC  (%) 

LogRandom  (%) 

LHS  (%) 

CVT  (%) 

Blade  1  WPC  mean 

0.00789 

0.01298 

0.0231 

0.0176 

Blade  1  WPC  variance 

0.2 

1.44 

0.657 

1.544 

Blade  2  WPC  mean 

0.003 

0.01299 

0.02113 

0.0132 

Blade  2  WPC  variance 

0.083 

0.16 

0.647 

0.315 

Table  6.5:  Model-constrained  (MC)  adaptive  sampling  method  versus  logarithmic 
random  (LogRandom)  sampling,  LHS,  and  CVT  for  the  case  with  ten  geometric 
variability  parameters.  For  all  methods,  29  sample  points  are  generated,  and  the 
percentage  errors  in  predictions  of  means  and  variances  of  work  per  cycle  (WPC)  of 
the  resulting  reduced  models  are  shown  for  10,000  randomly  selected  blade  geometries. 


MC  (%) 

LogRandom  (%) 

LHS  (%) 

CVT  (%) 

Blade  1  WPC  mean 

0.3682 

0.0428 

0.025 

0.029 

Blade  1  WPC  variance 

0.6082 

3.4950 

0.3211 

0.3759 

Blade  2  WPC  mean 

0.0869 

0.091 

0.0054 

0.0209 

Blade  2  WPC  variance 

1.879 

4.819 

2.435 

0.055 
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Chapter  7 


Conclusions  and  Recommendations 

7.1  Thesis  Summary  and  Contributions 

In  this  thesis,  we  have  proposed  a  model-constrained  greedy-based  adaptive  sampling 
approach  to  address  the  challenge  of  sampling  high-dimensional  parameter  spaces  for 
model  reduction  of  large-scale  systems.  The  method  provides  a  systematic  procedure 
for  sampling  the  parametric  input  space,  and  the  forcing  input  space.  In  particular, 
we  formulate  the  sampling  problem  as  an  optimization  problem  that  targets  an  error 
estimation  (which  could  be  an  output  error  indicator,  an  output  error  bound,  or  the 
true  output  error)  of  reduced  model  output  prediction.  The  optimization  problem 
is  defined  by  introducing  as  constraints  the  systems  of  equations  representing  the 
reduced  model  (and  possibly  the  full  model  if  the  true  output  error  is  targeted).  The 
optimization  formulation  treats  the  parameter  space  as  continuous;  that  is,  we  do 
not  require  a  priori  selection  of  a  discrete  training  parameter  set.  Further,  since  any 
error  estimation  can  be  used  as  our  selection  criteria,  our  approach  is  applicable  in 
cases  for  which  output  error  bounds  are  unavailable.  Finally,  we  use  a  state-of-the- 
art  optimization  technique,  namely  the  subspace  trust-region  interior  reflective  inex¬ 
act  Newton-CG  method,  to  solve  the  resulting  greedy  PDE-constrained  optimization 
problem. 

In  principle,  one  can  treat  the  forcing  input  as  another  parameter,  albeit  a  special 
parameter,  and  the  adaptive  sampling  method  can  be  applied  for  both  parametric 
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input  space  and  unsteady  forcing  space. 

We  have  also  provided  an  analysis  of  the  model-constrained  adaptive  sampling 
approaches  with  two  key  results  as  follows.  First,  the  greedy  optimization  problem 
is  well-posed  in  the  sense  that  there  exists  at  least  one  solution.  Second,  if  the  full 
model  is  linear  in  the  state  vector,  then  the  adaptive  sampling  approach  will  never 
sample  at  the  previous  sampled  points  in  the  parameter  space. 

The  model-constrained  adaptive  sampling  approach  has  been  applied  to  a  steady 
thermal  fin  optimal  design  problem  and  to  probabilistic  analysis  of  geometric  variabil¬ 
ity  in  a  turbomachinery  application.  While  these  two  examples  are  both  linear  in  the 
state  vector,  our  sampling  approach  could  also  be  applied  to  general  nonlinear  prob¬ 
lems.  In  the  nonlinear  case,  efficiency  of  the  reduced  model  could  be  addressed,  for 
example,  using  the  interpolation  method  of  [61,62,68].  In  the  following,  we  summarize 
the  results  and  the  conclusions  from  numerical  studies  of  these  two  examples. 

The  model-constrained  adaptive  reduction  method  has  been  able  to  create  reduced 
models  that  accurately  represent  the  full  models  over  a  wide  range  of  parameter 
values  in  high-dimensional  parametric  spaces.  Since  the  thermal  fin  problem  is  quite 
small,  the  offline  cost  of  generating  the  reduced  model  is  not  offset  by  the  savings 
in  online  cost  even  though  using  the  reduced  model  in  optimal  design  is  an  order 
of  magnitude  faster  than  using  the  full  model.  The  trade-off  between  the  online  and 
offline  cost  has  been  seen  to  be  much  more  dramatic  for  probabilistic  analysis  problems 
in  which  the  savings  in  the  online  cost  are  substantial  and  offset  the  offline  cost. 
In  particular,  for  probabilistic  analysis  of  a  subsonic  blade  row  with  ten  geometric 
variability  parameters,  the  reduced  model  provided  a  factor  of  468  speed-up  for  a 
Monte  Carlo  simulation  of  10,000  samples. 

Either  the  true  error  and  the  error  indicator  can  be  used  as  sampling  selection 
in  the  model-constrained  adaptive  sampling.  Which  method  to  use  is  a  matter  of 
trading  offline  cost  with  reduced  model  accuracy.  For  the  thermal  fin  problem,  using 
the  error- indicator  approach  leads  to  less  expensive  offline  cost,  but  the  resulting 
reduced  model  can  be  less  accurate  than  that  obtained  using  the  true-error  approach 
with  the  same  basis  size.  However,  both  steady  and  unsteady  results  have  shown 
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that,  using  the  error  indicator,  in  particular  the  residual,  one  can  obtain  reduced 
models  that  are  accurate  over  a  wide  range  of  parameter  values  in  multi- dimensional 
parametric  input  spaces. 

Both  steady  and  unsteady  results  numerically  show  that  the  subspace  trust-region 
interior  reflective  inexact  Newton  CG  solver  is  efficient  for  the  greedy  optimization 
problem  in  the  sense  that  the  number  of  Newton  steps  scales  as  0(d ),  where  d  is 
the  number  of  parameters.  In  addition,  one  can  also  theoretically  show  that  if  a 
good  preconditioner  for  the  CG  solver  is  available  so  that  the  number  of  Hessian- 
vector  products  is  independent  of  the  number  of  the  parameters,  the  offline  cost  for 
each  greedy  cycle  scales  linearly  with  the  dimension  of  the  parameter  space,  d.  In 
addition,  if  the  number  of  greedy  cycles  is  fixed,  i.e.  the  size  of  the  reduced  model  is 
given,  then  one  can  prove  that  the  total  offline  cost  of  constructing  the  reduced  basis 
scales  linearly  with  the  number  of  parameters.  As  a  result,  the  model-constrained 
adaptive  sampling  approach  can  be  used  for  problems  that  depend  on  a  large  number 
of  parameters. 

In  this  thesis,  various  sampling  methods  are  compared  in  the  context  of  snapshot 
generation  for  model  reduction.  The  numerical  results  for  the  thermal  fin  show  that 
with  the  same  number  of  sampling  points,  the  resulting  reduced  model  from  logarith¬ 
mic  random  sampling  approach  is  more  accurate  than  those  obtained  from  LHS  and 
CVT  sampling  methods.  Compared  to  these  methods,  the  model-constrained  sam¬ 
pling  yields  a  reduced  model  with  error  of  several  orders  of  magnitude  smaller  than 
those  obtained  using  logarithmic  random  sampling,  uniform  random  sampling,  LHS 
and  CVT  sampling  methods.  As  the  size  of  the  parametric  input  space  increases, 
the  difference  in  accuracy  of  the  reduced  model  using  model-constrained  sampling 
method  and  those  of  other  methods  is  larger.  This  reflects  the  fact  that  the  model- 
constrained  sampling  is  a  model-based  sampling  method;  that  is,  the  parameter  is 
sampled  where  the  indicator  of  the  error  between  the  full  and  the  reduced  models  is 
locally  largest,  whereas  the  other  statistically-based  sampling  methods  use  no  knowl¬ 
edge  of  the  underlying  model  when  selecting  their  sample  points.  As  the  dimension  of 
the  parametric  space  increases,  it  becomes  more  and  more  difficult  to  adequately  cover 
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the  space  with  a  reasonable  number  of  samples  using  the  statistically-based  sampling 
methods.  However,  if  the  problems  under  consideration  are  weakly  nonlinear  in  the 
parameters  or  the  parameter  ranges  are  small,  the  results  of  the  linearized  CFD  model 
show  that  the  model-constrained  sampling  approach  and  the  other  statistically-based 
sampling  methods  yield  sample  points  of  similar  quality  in  the  context  of  snapshot 
generation  for  model  reduction. 


7.2  Extensions  and  Future  Work 

There  are  a  number  of  extensions  that  could  be  applied  to  the  model-constrained 
greedy-based  adaptive  sampling  approach  in  this  research.  These  extensions  include 
improvements  to  the  adaptive  sampling  method,  possible  treatments  for  its  limita¬ 
tions,  and  possibilities  for  additional  model  applications.  In  the  following,  we  present 
some  extensions  and  future  work. 

We  have  used  a  “black-box”  approach  (also  called  the  reduced-space  approach)  to 
solve  the  greedy  optimization  problem.  That  is,  the  states  are  eliminated  by  solving 
the  PDE  constraints  exactly  (up  to  machine  zero),  and  hence  the  PDE  constraints 
are  not  visible  to  the  optimizer.  As  a  result,  feasibility  is  always  guaranteed  and 
the  black-box  approach  only  moves  towards  optimality.  The  all-at-once  approach 
(sometimes  called  the  full-space  approach)  is  more  efficient  in  the  sense  that  it  moves 
towards  optimality  and  feasibility  at  the  same  time;  hence  it  does  not  require  the 
possibly  very  expensive  solution  of  the  PDE  constraints,  especially  for  problems  that 
are  nonlinear  in  the  state  vector. 

It  should  be  emphasized  again  that  while  we  have  mainly  presented  both  the¬ 
oretical  and  numerical  results  for  problems  that  are  linear  in  the  state  vector,  the 
model-constrained  sampling  approach  is  applicable  for  problems  that  are  nonlinear  in 
the  state  vector  as  well.  Therefore,  if  a  numerical  solver  is  available  for  the  nonlinear 
problem  under  consideration — in  particular,  given  a  parameter  set  as  the  input,  the 
output  of  the  solver  is  the  state  solution  corresponding  to  that  input  set — it  is  straight¬ 
forward  to  apply  the  adaptive  sampling  approach  using  the  framework  proposed  in 
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this  thesis.  However,  the  main  assumption  in  Corollary  3.2  is  that  the  problem  under 
consideration  is  linear  in  the  state  vector  so  that  the  a  priori  result  of  the  form  stated 
in  Theorem  2.2  holds.  Note  that  if  a  priori  result  of  this  form  is  also  available  for 
a  problem  that  is  nonlinear  in  the  state  vector,  Corollary  3.2  still  holds.  This  is, 
however,  not  a  severe  limitation  if,  for  example,  one  uses  the  empirical  interpolation 
method  [61,62,68]  and  its  recent  extensions  [69,114,115]  to  pre-process  the  nonlin¬ 
ear  terms  by  a  linear  combination  of  empirical  interpolation  basis  functions.  In  that 
case,  the  results  in  Corollary  3.2  still  apply  for  the  pre-processed  model.  It  should 
be  pointed  out  that  while  the  model-constrained  sampling  approach  addresses  the 
question  on  how  to  construct  a  reduced  basis  that  is  valid  over  wide  ranges  of  param¬ 
eter  values  in  multi- dimensional  parametric  input  space,  the  empirical  interpolation 
approach  addresses  the  question  on  how  to  construct  an  efficient  reduced  model  given 
a  reduced  basis.  Therefore,  the  combination  of  the  model-constrained  approach  and 
the  empirical  interpolation  method  is  a  promising  model  reduction  method  that  yields 
accurate  and  efficient  reduced  models  for  a  general  nonlinear  problem  that  depends  on 
a  large  number  of  parameters.  An  example  in  which  this  combination  could  be  used 
is  the  probabilistic  analysis.  For  probabilistic  analysis  applications  in  this  research, 
we  have  limited  ourselves  to  a  linearized  model  in  both  states  and  parameters  so  that 
efficient  reduced  models  are  available  for  the  online  stage.  Therefore  the  model  is 
only  valid  for  a  small  perturbation  in  states  and  parameters  about  the  linearization 
point.  An  improved  model  is  therefore  to  combine  the  adaptive  sampling  method  and 
the  empirical  interpolation  method  so  that  the  reduced  model  is  applicable  for  wider 
ranges  of  geometric  variations. 
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